|
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, class ArrayOutFields, class ArrayInFieldsLeft, class ArrayInFieldsRight> 00038 void ArrayTools::contractFieldFieldScalar(ArrayOutFields & outputFields, 00039 const ArrayInFieldsLeft & leftFields, 00040 const ArrayInFieldsRight & rightFields, 00041 const ECompEngine compEngine, 00042 const bool sumInto) { 00043 00044 00045 #ifdef HAVE_INTREPID_DEBUG 00046 TEST_FOR_EXCEPTION( (leftFields.rank() != 3 ), std::invalid_argument, 00047 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Rank of the left input argument must equal 3!"); 00048 TEST_FOR_EXCEPTION( (rightFields.rank() != 3 ), std::invalid_argument, 00049 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Rank of right input argument must equal 3!"); 00050 TEST_FOR_EXCEPTION( (outputFields.rank() != 3 ), std::invalid_argument, 00051 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Rank of output argument must equal 3!"); 00052 TEST_FOR_EXCEPTION( (leftFields.dimension(0) != rightFields.dimension(0) ), std::invalid_argument, 00053 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!"); 00054 TEST_FOR_EXCEPTION( (leftFields.dimension(2) != rightFields.dimension(2) ), std::invalid_argument, 00055 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Second dimensions (numbers of integration points) of the left and right input containers must agree!"); 00056 TEST_FOR_EXCEPTION( (outputFields.dimension(0) != rightFields.dimension(0) ), std::invalid_argument, 00057 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!"); 00058 TEST_FOR_EXCEPTION( (outputFields.dimension(1) != leftFields.dimension(1) ), std::invalid_argument, 00059 ">>> ERROR (ArrayTools::contractFieldFieldScalar): First dimension of output container and first dimension of left input container must agree!"); 00060 TEST_FOR_EXCEPTION( (outputFields.dimension(2) != rightFields.dimension(1) ), std::invalid_argument, 00061 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Second dimension of output container and first dimension of right input container must agree!"); 00062 #endif 00063 00064 // get sizes 00065 int numCells = leftFields.dimension(0); 00066 int numLeftFields = leftFields.dimension(1); 00067 int numRightFields = rightFields.dimension(1); 00068 int numPoints = leftFields.dimension(2); 00069 00070 switch(compEngine) { 00071 case COMP_CPP: { 00072 if (sumInto) { 00073 for (int cl = 0; cl < numCells; cl++) { 00074 for (int lbf = 0; lbf < numLeftFields; lbf++) { 00075 for (int rbf = 0; rbf < numRightFields; rbf++) { 00076 Scalar tmpVal(0); 00077 for (int qp = 0; qp < numPoints; qp++) { 00078 tmpVal += leftFields(cl, lbf, qp)*rightFields(cl, rbf, qp); 00079 } // P-loop 00080 outputFields(cl, lbf, rbf) += tmpVal; 00081 } // R-loop 00082 } // L-loop 00083 } // C-loop 00084 } 00085 else { 00086 for (int cl = 0; cl < numCells; cl++) { 00087 for (int lbf = 0; lbf < numLeftFields; lbf++) { 00088 for (int rbf = 0; rbf < numRightFields; rbf++) { 00089 Scalar tmpVal(0); 00090 for (int qp = 0; qp < numPoints; qp++) { 00091 tmpVal += leftFields(cl, lbf, qp)*rightFields(cl, rbf, qp); 00092 } // P-loop 00093 outputFields(cl, lbf, rbf) = tmpVal; 00094 } // R-loop 00095 } // L-loop 00096 } // C-loop 00097 } 00098 } 00099 break; 00100 00101 case COMP_BLAS: { 00102 /* 00103 GEMM parameters and their values. 00104 (Note: It is assumed that the result needs to be transposed into row-major format. 00105 Think of left and right input matrices as A(p x l) and B(p x r), respectively, 00106 even though the indexing is ((C),L,P) and ((C),R,P). Due to BLAS formatting 00107 assumptions, we are computing (A^T*B)^T = B^T*A.) 00108 TRANSA TRANS 00109 TRANSB NO_TRANS 00110 M #rows(B^T) = number of right fields 00111 N #cols(A) = number of left fields 00112 K #cols(B^T) = number of integration points 00113 ALPHA 1.0 00114 A right data for cell cl = &rightFields[cl*skipR] 00115 LDA #rows(B) = number of integration points 00116 B left data for cell cl = &leftFields[cl*skipL] 00117 LDB #rows(A) = number of integration points 00118 BETA 0.0 00119 C result for cell cl = &outputFields[cl*skipOp] 00120 LDC #rows(C) = number of right fields 00121 */ 00122 int skipL = numLeftFields*numPoints; // size of the left data chunk per cell 00123 int skipR = numRightFields*numPoints; // size of the right data chunk per cell 00124 int skipOp = numLeftFields*numRightFields; // size of the output data chunk per cell 00125 Scalar alpha(1.0); // these are left unchanged by GEMM 00126 Scalar beta(0.0); 00127 if (sumInto) { 00128 beta = 1.0; 00129 } 00130 00131 for (int cl=0; cl < numCells; cl++) { 00132 /* Use this if data is used in row-major format */ 00133 Teuchos::BLAS<int, Scalar> myblas; 00134 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS, 00135 numRightFields, numLeftFields, numPoints, 00136 alpha, &rightFields[cl*skipR], numPoints, 00137 &leftFields[cl*skipL], numPoints, 00138 beta, &outputFields[cl*skipOp], numRightFields); 00139 /* Use this if data is used in column-major format */ 00140 /* 00141 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS, 00142 numLeftFields, numRightFields, numPoints, 00143 alpha, &leftFields[cl*skipL], numPoints, 00144 &rightFields[cl*skipR], numPoints, 00145 beta, &outputFields[cl*skipOp], numLeftFields); 00146 */ 00147 } 00148 } 00149 break; 00150 00151 default: 00152 TEST_FOR_EXCEPTION( ( ~isValidCompEngine(compEngine) ), std::invalid_argument, 00153 ">>> ERROR (ArrayTools::contractFieldFieldScalar): Computational engine not defined!"); 00154 } // switch(compEngine) 00155 } // end contractFieldFieldScalar 00156 00157 00158 template<class Scalar, class ArrayOutFields, class ArrayInFieldsLeft, class ArrayInFieldsRight> 00159 void ArrayTools::contractFieldFieldVector(ArrayOutFields & outputFields, 00160 const ArrayInFieldsLeft & leftFields, 00161 const ArrayInFieldsRight & rightFields, 00162 const ECompEngine compEngine, 00163 const bool sumInto) { 00164 00165 #ifdef HAVE_INTREPID_DEBUG 00166 TEST_FOR_EXCEPTION( (leftFields.rank() != 4 ), std::invalid_argument, 00167 ">>> ERROR (ArrayTools::contractFieldFieldVector): Rank of the left input argument must equal 4!"); 00168 TEST_FOR_EXCEPTION( (rightFields.rank() != 4 ), std::invalid_argument, 00169 ">>> ERROR (ArrayTools::contractFieldFieldVector): Rank of right input argument must equal 4!"); 00170 TEST_FOR_EXCEPTION( (outputFields.rank() != 3 ), std::invalid_argument, 00171 ">>> ERROR (ArrayTools::contractFieldFieldVector): Rank of output argument must equal 3!"); 00172 TEST_FOR_EXCEPTION( (leftFields.dimension(0) != rightFields.dimension(0) ), std::invalid_argument, 00173 ">>> ERROR (ArrayTools::contractFieldFieldVector): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!"); 00174 TEST_FOR_EXCEPTION( (leftFields.dimension(2) != rightFields.dimension(2) ), std::invalid_argument, 00175 ">>> ERROR (ArrayTools::contractFieldFieldVector): Second dimensions (numbers of integration points) of the left and right input containers must agree!"); 00176 TEST_FOR_EXCEPTION( (leftFields.dimension(3) != rightFields.dimension(3) ), std::invalid_argument, 00177 ">>> ERROR (ArrayTools::contractFieldFieldVector): Third dimensions (numbers of vector components) of the left and right input containers must agree!"); 00178 TEST_FOR_EXCEPTION( (outputFields.dimension(0) != rightFields.dimension(0) ), std::invalid_argument, 00179 ">>> ERROR (ArrayTools::contractFieldFieldVector): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!"); 00180 TEST_FOR_EXCEPTION( (outputFields.dimension(1) != leftFields.dimension(1) ), std::invalid_argument, 00181 ">>> ERROR (ArrayTools::contractFieldFieldVector): First dimension of output container and first dimension of left input container must agree!"); 00182 TEST_FOR_EXCEPTION( (outputFields.dimension(2) != rightFields.dimension(1) ), std::invalid_argument, 00183 ">>> ERROR (ArrayTools::contractFieldFieldVector): Second dimension of output container and first dimension of right input container must agree!"); 00184 #endif 00185 00186 // get sizes 00187 int numCells = leftFields.dimension(0); 00188 int numLeftFields = leftFields.dimension(1); 00189 int numRightFields = rightFields.dimension(1); 00190 int numPoints = leftFields.dimension(2); 00191 int dimVec = leftFields.dimension(3); 00192 00193 switch(compEngine) { 00194 case COMP_CPP: { 00195 if (sumInto) { 00196 for (int cl = 0; cl < numCells; cl++) { 00197 for (int lbf = 0; lbf < numLeftFields; lbf++) { 00198 for (int rbf = 0; rbf < numRightFields; rbf++) { 00199 Scalar tmpVal(0); 00200 for (int qp = 0; qp < numPoints; qp++) { 00201 for (int iVec = 0; iVec < dimVec; iVec++) { 00202 tmpVal += leftFields(cl, lbf, qp, iVec)*rightFields(cl, rbf, qp, iVec); 00203 } //D-loop 00204 } // P-loop 00205 outputFields(cl, lbf, rbf) += tmpVal; 00206 } // R-loop 00207 } // L-loop 00208 } // C-loop 00209 } 00210 else { 00211 for (int cl = 0; cl < numCells; cl++) { 00212 for (int lbf = 0; lbf < numLeftFields; lbf++) { 00213 for (int rbf = 0; rbf < numRightFields; rbf++) { 00214 Scalar tmpVal(0); 00215 for (int qp = 0; qp < numPoints; qp++) { 00216 for (int iVec = 0; iVec < dimVec; iVec++) { 00217 tmpVal += leftFields(cl, lbf, qp, iVec)*rightFields(cl, rbf, qp, iVec); 00218 } //D-loop 00219 } // P-loop 00220 outputFields(cl, lbf, rbf) = tmpVal; 00221 } // R-loop 00222 } // L-loop 00223 } // C-loop 00224 } 00225 } 00226 break; 00227 00228 case COMP_BLAS: { 00229 /* 00230 GEMM parameters and their values. 00231 (Note: It is assumed that the result needs to be transposed into row-major format. 00232 Think of left and right input matrices as A(p x l) and B(p x r), respectively, 00233 even though the indexing is ((C),L,P) and ((C),R,P). Due to BLAS formatting 00234 assumptions, we are computing (A^T*B)^T = B^T*A.) 00235 TRANSA TRANS 00236 TRANSB NO_TRANS 00237 M #rows(B^T) = number of right fields 00238 N #cols(A) = number of left fields 00239 K #cols(B^T) = number of integration points * size of vector 00240 ALPHA 1.0 00241 A right data for cell cl = &rightFields[cl*skipR] 00242 LDA #rows(B) = number of integration points * size of vector 00243 B left data for cell cl = &leftFields[cl*skipL] 00244 LDB #rows(A) = number of integration points * size of vector 00245 BETA 0.0 00246 C result for cell cl = &outputFields[cl*skipOp] 00247 LDC #rows(C) = number of right fields 00248 */ 00249 int numData = numPoints*dimVec; 00250 int skipL = numLeftFields*numData; // size of the left data chunk per cell 00251 int skipR = numRightFields*numData; // size of the right data chunk per cell 00252 int skipOp = numLeftFields*numRightFields; // size of the output data chunk per cell 00253 Scalar alpha(1.0); // these are left unchanged by GEMM 00254 Scalar beta(0.0); 00255 if (sumInto) { 00256 beta = 1.0; 00257 } 00258 00259 for (int cl=0; cl < numCells; cl++) { 00260 /* Use this if data is used in row-major format */ 00261 Teuchos::BLAS<int, Scalar> myblas; 00262 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS, 00263 numRightFields, numLeftFields, numData, 00264 alpha, &rightFields[cl*skipR], numData, 00265 &leftFields[cl*skipL], numData, 00266 beta, &outputFields[cl*skipOp], numRightFields); 00267 /* Use this if data is used in column-major format */ 00268 /* 00269 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS, 00270 numLeftFields, numRightFields, numData, 00271 alpha, &leftFields[cl*skipL], numData, 00272 &rightFields[cl*skipR], numData, 00273 beta, &outputFields[cl*skipOp], numLeftFields); 00274 */ 00275 } 00276 } 00277 break; 00278 00279 default: 00280 TEST_FOR_EXCEPTION( ( ~isValidCompEngine(compEngine) ), std::invalid_argument, 00281 ">>> ERROR (ArrayTools::contractFieldFieldVector): Computational engine not defined!"); 00282 } // switch(compEngine) 00283 } // end contractFieldFieldVector 00284 00285 00286 template<class Scalar, class ArrayOutFields, class ArrayInFieldsLeft, class ArrayInFieldsRight> 00287 void ArrayTools::contractFieldFieldTensor(ArrayOutFields & outputFields, 00288 const ArrayInFieldsLeft & leftFields, 00289 const ArrayInFieldsRight & rightFields, 00290 const ECompEngine compEngine, 00291 const bool sumInto) { 00292 00293 #ifdef HAVE_INTREPID_DEBUG 00294 TEST_FOR_EXCEPTION( (leftFields.rank() != 5 ), std::invalid_argument, 00295 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Rank of the left input argument must equal 5!"); 00296 TEST_FOR_EXCEPTION( (rightFields.rank() != 5 ), std::invalid_argument, 00297 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Rank of right input argument must equal 5!"); 00298 TEST_FOR_EXCEPTION( (outputFields.rank() != 3 ), std::invalid_argument, 00299 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Rank of output argument must equal 3!"); 00300 TEST_FOR_EXCEPTION( (leftFields.dimension(0) != rightFields.dimension(0) ), std::invalid_argument, 00301 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!"); 00302 TEST_FOR_EXCEPTION( (leftFields.dimension(2) != rightFields.dimension(2) ), std::invalid_argument, 00303 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Second dimensions (numbers of integration points) of the left and right input containers must agree!"); 00304 TEST_FOR_EXCEPTION( (leftFields.dimension(3) != rightFields.dimension(3) ), std::invalid_argument, 00305 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Third dimensions (first tensor dimensions) of the left and right input containers must agree!"); 00306 TEST_FOR_EXCEPTION( (leftFields.dimension(4) != rightFields.dimension(4) ), std::invalid_argument, 00307 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Fourth dimensions (second tensor dimensions) of the left and right input containers must agree!"); 00308 TEST_FOR_EXCEPTION( (outputFields.dimension(0) != rightFields.dimension(0) ), std::invalid_argument, 00309 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!"); 00310 TEST_FOR_EXCEPTION( (outputFields.dimension(1) != leftFields.dimension(1) ), std::invalid_argument, 00311 ">>> ERROR (ArrayTools::contractFieldFieldTensor): First dimension of output container and first dimension of left input container must agree!"); 00312 TEST_FOR_EXCEPTION( (outputFields.dimension(2) != rightFields.dimension(1) ), std::invalid_argument, 00313 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Second dimension of output container and first dimension of right input container must agree!"); 00314 #endif 00315 00316 // get sizes 00317 int numCells = leftFields.dimension(0); 00318 int numLeftFields = leftFields.dimension(1); 00319 int numRightFields = rightFields.dimension(1); 00320 int numPoints = leftFields.dimension(2); 00321 int dim1Tensor = leftFields.dimension(3); 00322 int dim2Tensor = leftFields.dimension(4); 00323 00324 switch(compEngine) { 00325 case COMP_CPP: { 00326 if (sumInto) { 00327 for (int cl = 0; cl < numCells; cl++) { 00328 for (int lbf = 0; lbf < numLeftFields; lbf++) { 00329 for (int rbf = 0; rbf < numRightFields; rbf++) { 00330 Scalar tmpVal(0); 00331 for (int qp = 0; qp < numPoints; qp++) { 00332 for (int iTens1 = 0; iTens1 < dim1Tensor; iTens1++) { 00333 for (int iTens2 = 0; iTens2 < dim2Tensor; iTens2++) { 00334 tmpVal += leftFields(cl, lbf, qp, iTens1, iTens2)*rightFields(cl, rbf, qp, iTens1, iTens2); 00335 } // D2-loop 00336 } // D1-loop 00337 } // P-loop 00338 outputFields(cl, lbf, rbf) += tmpVal; 00339 } // R-loop 00340 } // L-loop 00341 } // C-loop 00342 } 00343 else { 00344 for (int cl = 0; cl < numCells; cl++) { 00345 for (int lbf = 0; lbf < numLeftFields; lbf++) { 00346 for (int rbf = 0; rbf < numRightFields; rbf++) { 00347 Scalar tmpVal(0); 00348 for (int qp = 0; qp < numPoints; qp++) { 00349 for (int iTens1 = 0; iTens1 < dim1Tensor; iTens1++) { 00350 for (int iTens2 = 0; iTens2 < dim2Tensor; iTens2++) { 00351 tmpVal += leftFields(cl, lbf, qp, iTens1, iTens2)*rightFields(cl, rbf, qp, iTens1, iTens2); 00352 } // D2-loop 00353 } // D1-loop 00354 } // P-loop 00355 outputFields(cl, lbf, rbf) = tmpVal; 00356 } // R-loop 00357 } // L-loop 00358 } // C-loop 00359 } 00360 } 00361 break; 00362 00363 case COMP_BLAS: { 00364 /* 00365 GEMM parameters and their values. 00366 (Note: It is assumed that the result needs to be transposed into row-major format. 00367 Think of left and right input matrices as A(p x l) and B(p x r), respectively, 00368 even though the indexing is ((C),L,P) and ((C),R,P). Due to BLAS formatting 00369 assumptions, we are computing (A^T*B)^T = B^T*A.) 00370 TRANSA TRANS 00371 TRANSB NO_TRANS 00372 M #rows(B^T) = number of right fields 00373 N #cols(A) = number of left fields 00374 K #cols(B^T) = number of integration points * size of tensor 00375 ALPHA 1.0 00376 A right data for cell cl = &rightFields[cl*skipR] 00377 LDA #rows(B) = number of integration points * size of tensor 00378 B left data for cell cl = &leftFields[cl*skipL] 00379 LDB #rows(A) = number of integration points * size of tensor 00380 BETA 0.0 00381 C result for cell cl = &outputFields[cl*skipOp] 00382 LDC #rows(C) = number of right fields 00383 */ 00384 int numData = numPoints*dim1Tensor*dim2Tensor; 00385 int skipL = numLeftFields*numData; // size of the left data chunk per cell 00386 int skipR = numRightFields*numData; // size of the right data chunk per cell 00387 int skipOp = numLeftFields*numRightFields; // size of the output data chunk per cell 00388 Scalar alpha(1.0); // these are left unchanged by GEMM 00389 Scalar beta(0.0); 00390 if (sumInto) { 00391 beta = 1.0; 00392 } 00393 00394 for (int cl=0; cl < numCells; cl++) { 00395 /* Use this if data is used in row-major format */ 00396 Teuchos::BLAS<int, Scalar> myblas; 00397 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS, 00398 numRightFields, numLeftFields, numData, 00399 alpha, &rightFields[cl*skipR], numData, 00400 &leftFields[cl*skipL], numData, 00401 beta, &outputFields[cl*skipOp], numRightFields); 00402 /* Use this if data is used in column-major format */ 00403 /* 00404 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS, 00405 numLeftFields, numRightFields, numData, 00406 alpha, &leftFields[cl*skipL], numData, 00407 &rightFields[cl*skipR], numData, 00408 beta, &outputFields[cl*skipOp], numLeftFields); 00409 */ 00410 } 00411 } 00412 break; 00413 00414 default: 00415 TEST_FOR_EXCEPTION( ( ~isValidCompEngine(compEngine) ), std::invalid_argument, 00416 ">>> ERROR (ArrayTools::contractFieldFieldTensor): Computational engine not defined!"); 00417 } // switch(compEngine) 00418 } // end contractFieldFieldTensor 00419 00420 00421 template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields> 00422 void ArrayTools::contractDataFieldScalar(ArrayOutFields & outputFields, 00423 const ArrayInData & inputData, 00424 const ArrayInFields & inputFields, 00425 const ECompEngine compEngine, 00426 const bool sumInto) { 00427 00428 00429 #ifdef HAVE_INTREPID_DEBUG 00430 TEST_FOR_EXCEPTION( (inputFields.rank() != 3 ), std::invalid_argument, 00431 ">>> ERROR (ArrayTools::contractDataFieldScalar): Rank of the fields input argument must equal 3!"); 00432 TEST_FOR_EXCEPTION( (inputData.rank() != 2 ), std::invalid_argument, 00433 ">>> ERROR (ArrayTools::contractDataFieldScalar): Rank of the data input argument must equal 2!"); 00434 TEST_FOR_EXCEPTION( (outputFields.rank() != 2 ), std::invalid_argument, 00435 ">>> ERROR (ArrayTools::contractDataFieldScalar): Rank of output argument must equal 2!"); 00436 TEST_FOR_EXCEPTION( (inputFields.dimension(0) != inputData.dimension(0) ), std::invalid_argument, 00437 ">>> ERROR (ArrayTools::contractDataFieldScalar): Zeroth dimensions (number of integration domains) of the fields and data input containers must agree!"); 00438 TEST_FOR_EXCEPTION( ( (inputFields.dimension(2) != inputData.dimension(1)) && (inputData.dimension(1) != 1) ), std::invalid_argument, 00439 ">>> ERROR (ArrayTools::contractDataFieldScalar): Second dimension of fields input container and first dimension of data input container (number of integration points) must agree or first data dimension must be 1!"); 00440 TEST_FOR_EXCEPTION( (outputFields.dimension(0) != inputFields.dimension(0) ), std::invalid_argument, 00441 ">>> ERROR (ArrayTools::contractDataFieldScalar): Zeroth dimensions (numbers of integration domains) of the fields input and output containers must agree!"); 00442 TEST_FOR_EXCEPTION( (outputFields.dimension(1) != inputFields.dimension(1) ), std::invalid_argument, 00443 ">>> ERROR (ArrayTools::contractDataFieldScalar): First dimensions (number of fields) of the fields input and output containers must agree!"); 00444 #endif 00445 00446 // get sizes 00447 int numCells = inputFields.dimension(0); 00448 int numFields = inputFields.dimension(1); 00449 int numPoints = inputFields.dimension(2); 00450 int numDataPoints = inputData.dimension(1); 00451 00452 ECompEngine myCompEngine = (numDataPoints == 1 ? COMP_CPP : compEngine); 00453 00454 switch(myCompEngine) { 00455 case COMP_CPP: { 00456 if (sumInto) { 00457 if (numDataPoints != 1) { // nonconstant data 00458 for (int cl = 0; cl < numCells; cl++) { 00459 for (int lbf = 0; lbf < numFields; lbf++) { 00460 Scalar tmpVal(0); 00461 for (int qp = 0; qp < numPoints; qp++) { 00462 tmpVal += inputFields(cl, lbf, qp)*inputData(cl, qp); 00463 } // P-loop 00464 outputFields(cl, lbf) += tmpVal; 00465 } // F-loop 00466 } // C-loop 00467 } 00468 else { // constant data 00469 for (int cl = 0; cl < numCells; cl++) { 00470 for (int lbf = 0; lbf < numFields; lbf++) { 00471 Scalar tmpVal(0); 00472 for (int qp = 0; qp < numPoints; qp++) { 00473 tmpVal += inputFields(cl, lbf, qp)*inputData(cl, 0); 00474 } // P-loop 00475 outputFields(cl, lbf) += tmpVal; 00476 } // F-loop 00477 } // C-loop 00478 } // numDataPoints 00479 } 00480 else { 00481 if (numDataPoints != 1) { // nonconstant data 00482 for (int cl = 0; cl < numCells; cl++) { 00483 for (int lbf = 0; lbf < numFields; lbf++) { 00484 Scalar tmpVal(0); 00485 for (int qp = 0; qp < numPoints; qp++) { 00486 tmpVal += inputFields(cl, lbf, qp)*inputData(cl, qp); 00487 } // P-loop 00488 outputFields(cl, lbf) = tmpVal; 00489 } // F-loop 00490 } // C-loop 00491 } 00492 else { // constant data 00493 for (int cl = 0; cl < numCells; cl++) { 00494 for (int lbf = 0; lbf < numFields; lbf++) { 00495 Scalar tmpVal(0); 00496 for (int qp = 0; qp < numPoints; qp++) { 00497 tmpVal += inputFields(cl, lbf, qp)*inputData(cl, 0); 00498 } // P-loop 00499 outputFields(cl, lbf) = tmpVal; 00500 } // F-loop 00501 } // C-loop 00502 } // numDataPoints 00503 } 00504 } 00505 break; 00506 00507 case COMP_BLAS: { 00508 /* 00509 GEMM parameters and their values. 00510 (Note: It is assumed that the result needs to be transposed into row-major format. 00511 Think of left and right input matrices as A(p x f) and B(p x 1), respectively, 00512 even though the indexing is ((C),F,P) and ((C),P). Due to BLAS formatting 00513 assumptions, we are computing (A^T*B)^T = B^T*A.) 00514 TRANSA TRANS 00515 TRANSB NO_TRANS 00516 M #rows(B^T) = 1 00517 N #cols(A) = number of input fields 00518 K #cols(B^T) = number of integration points * size of data 00519 ALPHA 1.0 00520 A right data for cell cl = &rightFields[cl*skipR] 00521 LDA #rows(B) = number of integration points * size of data 00522 B left data for cell cl = &leftFields[cl*skipL] 00523 LDB #rows(A) = number of integration points * size of data 00524 BETA 0.0 00525 C result for cell cl = &outputFields[cl*skipOp] 00526 LDC #rows(C) = 1 00527 */ 00528 int numData = numPoints; 00529 int skipL = numFields*numPoints; // size of the left data chunk per cell 00530 int skipR = numPoints; // size of the right data chunk per cell 00531 int skipOp = numFields; // size of the output data chunk per cell 00532 Scalar alpha(1.0); // these are left unchanged by GEMM 00533 Scalar beta(0.0); 00534 if (sumInto) { 00535 beta = 1.0; 00536 } 00537 00538 for (int cl=0; cl < numCells; cl++) { 00539 /* Use this if data is used in row-major format */ 00540 Teuchos::BLAS<int, Scalar> myblas; 00541 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS, 00542 1, numFields, numData, 00543 alpha, &inputData[cl*skipR], numData, 00544 &inputFields[cl*skipL], numData, 00545 beta, &outputFields[cl*skipOp], 1); 00546 /* Use this if data is used in column-major format */ 00547 /* 00548 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS, 00549 numFields, 1, numData, 00550 alpha, &inputFields[cl*skipL], numData, 00551 &inputData[cl*skipR], numData, 00552 beta, &outputFields[cl*skipOp], numFields); 00553 */ 00554 } 00555 } 00556 break; 00557 00558 default: 00559 TEST_FOR_EXCEPTION( ( ~isValidCompEngine(compEngine) ), std::invalid_argument, 00560 ">>> ERROR (ArrayTools::contractDataFieldScalar): Computational engine not defined!"); 00561 } // switch(compEngine) 00562 } // end contractDataFieldScalar 00563 00564 00565 template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields> 00566 void ArrayTools::contractDataFieldVector(ArrayOutFields & outputFields, 00567 const ArrayInData & inputData, 00568 const ArrayInFields & inputFields, 00569 const ECompEngine compEngine, 00570 const bool sumInto) { 00571 00572 #ifdef HAVE_INTREPID_DEBUG 00573 TEST_FOR_EXCEPTION( (inputFields.rank() != 4 ), std::invalid_argument, 00574 ">>> ERROR (ArrayTools::contractDataFieldVector): Rank of the fields input argument must equal 4!"); 00575 TEST_FOR_EXCEPTION( (inputData.rank() != 3 ), std::invalid_argument, 00576 ">>> ERROR (ArrayTools::contractDataFieldVector): Rank of the data input argument must equal 3!"); 00577 TEST_FOR_EXCEPTION( (outputFields.rank() != 2 ), std::invalid_argument, 00578 ">>> ERROR (ArrayTools::contractDataFieldVector): Rank of output argument must equal 2!"); 00579 TEST_FOR_EXCEPTION( (inputFields.dimension(0) != inputData.dimension(0) ), std::invalid_argument, 00580 ">>> ERROR (ArrayTools::contractDataFieldVector): Zeroth dimensions (number of integration domains) of the fields and data input containers must agree!"); 00581 TEST_FOR_EXCEPTION( ( (inputFields.dimension(2) != inputData.dimension(1)) && (inputData.dimension(1) != 1) ), std::invalid_argument, 00582 ">>> ERROR (ArrayTools::contractDataFieldVector): Second dimension of the fields input container and first dimension of data input container (number of integration points) must agree or first data dimension must be 1!"); 00583 TEST_FOR_EXCEPTION( (inputFields.dimension(3) != inputData.dimension(2) ), std::invalid_argument, 00584 ">>> ERROR (ArrayTools::contractDataFieldVector): Third dimension of the fields input container and second dimension of data input container (vector index) must agree!"); 00585 TEST_FOR_EXCEPTION( (outputFields.dimension(0) != inputFields.dimension(0) ), std::invalid_argument, 00586 ">>> ERROR (ArrayTools::contractDataFieldVector): Zeroth dimensions (numbers of integration domains) of the fields input and output containers must agree!"); 00587 TEST_FOR_EXCEPTION( (outputFields.dimension(1) != inputFields.dimension(1) ), std::invalid_argument, 00588 ">>> ERROR (ArrayTools::contractDataFieldVector): First dimensions of output container and fields input container (number of fields) must agree!"); 00589 #endif 00590 00591 // get sizes 00592 int numCells = inputFields.dimension(0); 00593 int numFields = inputFields.dimension(1); 00594 int numPoints = inputFields.dimension(2); 00595 int dimVec = inputFields.dimension(3); 00596 int numDataPoints = inputData.dimension(1); 00597 00598 ECompEngine myCompEngine = (numDataPoints == 1 ? COMP_CPP : compEngine); 00599 00600 switch(myCompEngine) { 00601 case COMP_CPP: { 00602 if (sumInto) { 00603 if (numDataPoints != 1) { // nonconstant data 00604 for (int cl = 0; cl < numCells; cl++) { 00605 for (int lbf = 0; lbf < numFields; lbf++) { 00606 Scalar tmpVal(0); 00607 for (int qp = 0; qp < numPoints; qp++) { 00608 for (int iVec = 0; iVec < dimVec; iVec++) { 00609 tmpVal += inputFields(cl, lbf, qp, iVec)*inputData(cl, qp, iVec); 00610 } // D-loop 00611 } // P-loop 00612 outputFields(cl, lbf) += tmpVal; 00613 } // F-loop 00614 } // C-loop 00615 } 00616 else { // constant data 00617 for (int cl = 0; cl < numCells; cl++) { 00618 for (int lbf = 0; lbf < numFields; lbf++) { 00619 Scalar tmpVal(0); 00620 for (int qp = 0; qp < numPoints; qp++) { 00621 for (int iVec = 0; iVec < dimVec; iVec++) { 00622 tmpVal += inputFields(cl, lbf, qp, iVec)*inputData(cl, 0, iVec); 00623 } //D-loop 00624 } // P-loop 00625 outputFields(cl, lbf) += tmpVal; 00626 } // F-loop 00627 } // C-loop 00628 } // numDataPoints 00629 } 00630 else { 00631 if (numDataPoints != 1) { // nonconstant data 00632 for (int cl = 0; cl < numCells; cl++) { 00633 for (int lbf = 0; lbf < numFields; lbf++) { 00634 Scalar tmpVal(0); 00635 for (int qp = 0; qp < numPoints; qp++) { 00636 for (int iVec = 0; iVec < dimVec; iVec++) { 00637 tmpVal += inputFields(cl, lbf, qp, iVec)*inputData(cl, qp, iVec); 00638 } // D-loop 00639 } // P-loop 00640 outputFields(cl, lbf) = tmpVal; 00641 } // F-loop 00642 } // C-loop 00643 } 00644 else { // constant data 00645 for (int cl = 0; cl < numCells; cl++) { 00646 for (int lbf = 0; lbf < numFields; lbf++) { 00647 Scalar tmpVal(0); 00648 for (int qp = 0; qp < numPoints; qp++) { 00649 for (int iVec = 0; iVec < dimVec; iVec++) { 00650 tmpVal += inputFields(cl, lbf, qp, iVec)*inputData(cl, 0, iVec); 00651 } //D-loop 00652 } // P-loop 00653 outputFields(cl, lbf) = tmpVal; 00654 } // F-loop 00655 } // C-loop 00656 } // numDataPoints 00657 } 00658 } 00659 break; 00660 00661 case COMP_BLAS: { 00662 /* 00663 GEMM parameters and their values. 00664 (Note: It is assumed that the result needs to be transposed into row-major format. 00665 Think of left and right input matrices as A(p x f) and B(p x 1), respectively, 00666 even though the indexing is ((C),F,P) and ((C),P). Due to BLAS formatting 00667 assumptions, we are computing (A^T*B)^T = B^T*A.) 00668 TRANSA TRANS 00669 TRANSB NO_TRANS 00670 M #rows(B^T) = 1 00671 N #cols(A) = number of input fields 00672 K #cols(B^T) = number of integration points * size of data 00673 ALPHA 1.0 00674 A right data for cell cl = &rightFields[cl*skipR] 00675 LDA #rows(B) = number of integration points * size of data 00676 B left data for cell cl = &leftFields[cl*skipL] 00677 LDB #rows(A) = number of integration points * size of data 00678 BETA 0.0 00679 C result for cell cl = &outputFields[cl*skipOp] 00680 LDC #rows(C) = 1 00681 */ 00682 int numData = numPoints*dimVec; 00683 int skipL = numFields*numData; // size of the left data chunk per cell 00684 int skipR = numData; // size of the right data chunk per cell 00685 int skipOp = numFields; // size of the output data chunk per cell 00686 Scalar alpha(1.0); // these are left unchanged by GEMM 00687 Scalar beta(0.0); 00688 if (sumInto) { 00689 beta = 1.0; 00690 } 00691 00692 for (int cl=0; cl < numCells; cl++) { 00693 /* Use this if data is used in row-major format */ 00694 Teuchos::BLAS<int, Scalar> myblas; 00695 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS, 00696 1, numFields, numData, 00697 alpha, &inputData[cl*skipR], numData, 00698 &inputFields[cl*skipL], numData, 00699 beta, &outputFields[cl*skipOp], 1); 00700 /* Use this if data is used in column-major format */ 00701 /* 00702 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS, 00703 numFields, 1, numData, 00704 alpha, &inputFields[cl*skipL], numData, 00705 &inputData[cl*skipR], numData, 00706 beta, &outputFields[cl*skipOp], numFields); 00707 */ 00708 } 00709 } 00710 break; 00711 00712 default: 00713 TEST_FOR_EXCEPTION( ( ~isValidCompEngine(compEngine) ), std::invalid_argument, 00714 ">>> ERROR (ArrayTools::contractDataFieldVector): Computational engine not defined!"); 00715 } // switch(compEngine) 00716 } // end contractDataFieldVector 00717 00718 00719 template<class Scalar, class ArrayOutFields, class ArrayInData, class ArrayInFields> 00720 void ArrayTools::contractDataFieldTensor(ArrayOutFields & outputFields, 00721 const ArrayInData & inputData, 00722 const ArrayInFields & inputFields, 00723 const ECompEngine compEngine, 00724 const bool sumInto) { 00725 00726 #ifdef HAVE_INTREPID_DEBUG 00727 TEST_FOR_EXCEPTION( (inputFields.rank() != 5 ), std::invalid_argument, 00728 ">>> ERROR (ArrayTools::contractDataFieldTensor): Rank of the fields input argument must equal 5!"); 00729 TEST_FOR_EXCEPTION( (inputData.rank() != 4 ), std::invalid_argument, 00730 ">>> ERROR (ArrayTools::contractDataFieldTensor): Rank of the data input argument must equal 4!"); 00731 TEST_FOR_EXCEPTION( (outputFields.rank() != 2 ), std::invalid_argument, 00732 ">>> ERROR (ArrayTools::contractDataFieldTensor): Rank of output argument must equal 2!"); 00733 TEST_FOR_EXCEPTION( (inputFields.dimension(0) != inputData.dimension(0) ), std::invalid_argument, 00734 ">>> ERROR (ArrayTools::contractDataFieldTensor): Zeroth dimensions (number of integration domains) of the fields and data input containers must agree!"); 00735 TEST_FOR_EXCEPTION( ( (inputFields.dimension(2) != inputData.dimension(1)) && (inputData.dimension(1) != 1) ), std::invalid_argument, 00736 ">>> ERROR (ArrayTools::contractDataFieldTensor): Second dimension of the fields input container and first dimension of data input container (number of integration points) must agree or first data dimension must be 1!"); 00737 TEST_FOR_EXCEPTION( (inputFields.dimension(3) != inputData.dimension(2) ), std::invalid_argument, 00738 ">>> ERROR (ArrayTools::contractDataFieldTensor): Third dimension of the fields input container and second dimension of data input container (first tensor dimension) must agree!"); 00739 TEST_FOR_EXCEPTION( (inputFields.dimension(4) != inputData.dimension(3) ), std::invalid_argument, 00740 ">>> ERROR (ArrayTools::contractDataFieldTensor): Fourth dimension of the fields input container and third dimension of data input container (second tensor dimension) must agree!"); 00741 TEST_FOR_EXCEPTION( (outputFields.dimension(0) != inputFields.dimension(0) ), std::invalid_argument, 00742 ">>> ERROR (ArrayTools::contractDataFieldTensor): Zeroth dimensions (numbers of integration domains) of the fields input and output containers must agree!"); 00743 TEST_FOR_EXCEPTION( (outputFields.dimension(1) != inputFields.dimension(1) ), std::invalid_argument, 00744 ">>> ERROR (ArrayTools::contractDataFieldTensor): First dimensions (number of fields) of output container and fields input container must agree!"); 00745 #endif 00746 00747 // get sizes 00748 int numCells = inputFields.dimension(0); 00749 int numFields = inputFields.dimension(1); 00750 int numPoints = inputFields.dimension(2); 00751 int dim1Tens = inputFields.dimension(3); 00752 int dim2Tens = inputFields.dimension(4); 00753 int numDataPoints = inputData.dimension(1); 00754 00755 ECompEngine myCompEngine = (numDataPoints == 1 ? COMP_CPP : compEngine); 00756 00757 switch(myCompEngine) { 00758 case COMP_CPP: { 00759 if (sumInto) { 00760 if (numDataPoints != 1) { // nonconstant data 00761 for (int cl = 0; cl < numCells; cl++) { 00762 for (int lbf = 0; lbf < numFields; lbf++) { 00763 Scalar tmpVal(0); 00764 for (int qp = 0; qp < numPoints; qp++) { 00765 for (int iTens1 = 0; iTens1 < dim1Tens; iTens1++) { 00766 for (int iTens2 =0; iTens2 < dim2Tens; iTens2++) { 00767 tmpVal += inputFields(cl, lbf, qp, iTens1, iTens2)*inputData(cl, qp, iTens1, iTens2); 00768 } // D2-loop 00769 } // D1-loop 00770 } // P-loop 00771 outputFields(cl, lbf) += tmpVal; 00772 } // F-loop 00773 } // C-loop 00774 } 00775 else { // constant data 00776 for (int cl = 0; cl < numCells; cl++) { 00777 for (int lbf = 0; lbf < numFields; lbf++) { 00778 Scalar tmpVal(0); 00779 for (int qp = 0; qp < numPoints; qp++) { 00780 for (int iTens1 = 0; iTens1 < dim1Tens; iTens1++) { 00781 for (int iTens2 = 0; iTens2 < dim2Tens; iTens2++) { 00782 tmpVal += inputFields(cl, lbf, qp, iTens1, iTens2)*inputData(cl, 0, iTens1, iTens2); 00783 } // D2-loop 00784 } // D1-loop 00785 } // P-loop 00786 outputFields(cl, lbf) += tmpVal; 00787 } // F-loop 00788 } // C-loop 00789 } // numDataPoints 00790 } 00791 else { 00792 if (numDataPoints != 1) { // nonconstant data 00793 for (int cl = 0; cl < numCells; cl++) { 00794 for (int lbf = 0; lbf < numFields; lbf++) { 00795 Scalar tmpVal(0); 00796 for (int qp = 0; qp < numPoints; qp++) { 00797 for (int iTens1 = 0; iTens1 < dim1Tens; iTens1++) { 00798 for (int iTens2 =0; iTens2 < dim2Tens; iTens2++) { 00799 tmpVal += inputFields(cl, lbf, qp, iTens1, iTens2)*inputData(cl, qp, iTens1, iTens2); 00800 } // D2-loop 00801 } // D1-loop 00802 } // P-loop 00803 outputFields(cl, lbf) = tmpVal; 00804 } // F-loop 00805 } // C-loop 00806 } 00807 else { // constant data 00808 for (int cl = 0; cl < numCells; cl++) { 00809 for (int lbf = 0; lbf < numFields; lbf++) { 00810 Scalar tmpVal(0); 00811 for (int qp = 0; qp < numPoints; qp++) { 00812 for (int iTens1 = 0; iTens1 < dim1Tens; iTens1++) { 00813 for (int iTens2 = 0; iTens2 < dim2Tens; iTens2++) { 00814 tmpVal += inputFields(cl, lbf, qp, iTens1, iTens2)*inputData(cl, 0, iTens1, iTens2); 00815 } // D2-loop 00816 } // D1-loop 00817 } // P-loop 00818 outputFields(cl, lbf) = tmpVal; 00819 } // F-loop 00820 } // C-loop 00821 } // numDataPoints 00822 } 00823 } 00824 break; 00825 00826 case COMP_BLAS: { 00827 /* 00828 GEMM parameters and their values. 00829 (Note: It is assumed that the result needs to be transposed into row-major format. 00830 Think of left and right input matrices as A(p x f) and B(p x 1), respectively, 00831 even though the indexing is ((C),F,P) and ((C),P). Due to BLAS formatting 00832 assumptions, we are computing (A^T*B)^T = B^T*A.) 00833 TRANSA TRANS 00834 TRANSB NO_TRANS 00835 M #rows(B^T) = 1 00836 N #cols(A) = number of input fields 00837 K #cols(B^T) = number of integration points * size of data 00838 ALPHA 1.0 00839 A right data for cell cl = &rightFields[cl*skipR] 00840 LDA #rows(B) = number of integration points * size of data 00841 B left data for cell cl = &leftFields[cl*skipL] 00842 LDB #rows(A) = number of integration points * size of data 00843 BETA 0.0 00844 C result for cell cl = &outputFields[cl*skipOp] 00845 LDC #rows(C) = 1 00846 */ 00847 int numData = numPoints*dim1Tens*dim2Tens; 00848 int skipL = numFields*numData; // size of the left data chunk per cell 00849 int skipR = numData; // size of the right data chunk per cell 00850 int skipOp = numFields; // size of the output data chunk per cell 00851 Scalar alpha(1.0); // these are left unchanged by GEMM 00852 Scalar beta(0.0); 00853 if (sumInto) { 00854 beta = 1.0; 00855 } 00856 00857 for (int cl=0; cl < numCells; cl++) { 00858 /* Use this if data is used in row-major format */ 00859 Teuchos::BLAS<int, Scalar> myblas; 00860 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS, 00861 1, numFields, numData, 00862 alpha, &inputData[cl*skipR], numData, 00863 &inputFields[cl*skipL], numData, 00864 beta, &outputFields[cl*skipOp], 1); 00865 /* Use this if data is used in column-major format */ 00866 /* 00867 myblas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS, 00868 numFields, 1, numData, 00869 alpha, &inputFields[cl*skipL], numData, 00870 &inputData[cl*skipR], numData, 00871 beta, &outputFields[cl*skipOp], numFields); 00872 */ 00873 } 00874 } 00875 break; 00876 00877 default: 00878 TEST_FOR_EXCEPTION( ( ~isValidCompEngine(compEngine) ), std::invalid_argument, 00879 ">>> ERROR (ArrayTools::contractDataFieldTensor): Computational engine not defined!"); 00880 } // switch(compEngine) 00881 } // end contractDataFieldTensor 00882 00883 00884 template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight> 00885 void ArrayTools::contractDataDataScalar(ArrayOutData & outputData, 00886 const ArrayInDataLeft & inputDataLeft, 00887 const ArrayInDataRight & inputDataRight, 00888 const ECompEngine compEngine, 00889 const bool sumInto) { 00890 00891 #ifdef HAVE_INTREPID_DEBUG 00892 TEST_FOR_EXCEPTION( (inputDataLeft.rank() != 2 ), std::invalid_argument, 00893 ">>> ERROR (ArrayTools::contractDataDataScalar): Rank of the left input argument must equal 2!"); 00894 TEST_FOR_EXCEPTION( (inputDataRight.rank() != 2 ), std::invalid_argument, 00895 ">>> ERROR (ArrayTools::contractDataDataScalar): Rank of right input argument must equal 2!"); 00896 TEST_FOR_EXCEPTION( (outputData.rank() != 1 ), std::invalid_argument, 00897 ">>> ERROR (ArrayTools::contractDataDataScalar): Rank of output argument must equal 1!"); 00898 TEST_FOR_EXCEPTION( (inputDataLeft.dimension(0) != inputDataRight.dimension(0) ), std::invalid_argument, 00899 ">>> ERROR (ArrayTools::contractDataDataScalar): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!"); 00900 TEST_FOR_EXCEPTION( (inputDataLeft.dimension(1) != inputDataRight.dimension(1) ), std::invalid_argument, 00901 ">>> ERROR (ArrayTools::contractDataDataScalar): First dimensions (numbers of integration points) of the left and right input containers must agree!"); 00902 TEST_FOR_EXCEPTION( (outputData.dimension(0) != inputDataRight.dimension(0) ), std::invalid_argument, 00903 ">>> ERROR (ArrayTools::contractDataDataScalar): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!"); 00904 #endif 00905 00906 // get sizes 00907 int numCells = inputDataLeft.dimension(0); 00908 int numPoints = inputDataLeft.dimension(1); 00909 00910 switch(compEngine) { 00911 case COMP_CPP: { 00912 if (sumInto) { 00913 for (int cl = 0; cl < numCells; cl++) { 00914 Scalar tmpVal(0); 00915 for (int qp = 0; qp < numPoints; qp++) { 00916 tmpVal += inputDataLeft(cl, qp)*inputDataRight(cl, qp); 00917 } // P-loop 00918 outputData(cl) += tmpVal; 00919 } // C-loop 00920 } 00921 else { 00922 for (int cl = 0; cl < numCells; cl++) { 00923 Scalar tmpVal(0); 00924 for (int qp = 0; qp < numPoints; qp++) { 00925 tmpVal += inputDataLeft(cl, qp)*inputDataRight(cl, qp); 00926 } // P-loop 00927 outputData(cl) = tmpVal; 00928 } // C-loop 00929 } 00930 } 00931 break; 00932 00933 case COMP_BLAS: { 00934 int incr = 1; // increment 00935 if (sumInto) { 00936 for (int cl=0; cl < numCells; cl++) { 00937 Teuchos::BLAS<int, Scalar> myblas; 00938 outputData(cl) += myblas.DOT(numPoints, &inputDataLeft[cl*numPoints], incr, &inputDataRight[cl*numPoints], incr); 00939 } 00940 } 00941 else { 00942 for (int cl=0; cl < numCells; cl++) { 00943 Teuchos::BLAS<int, Scalar> myblas; 00944 outputData(cl) = myblas.DOT(numPoints, &inputDataLeft[cl*numPoints], incr, &inputDataRight[cl*numPoints], incr); 00945 } 00946 } 00947 } 00948 break; 00949 00950 default: 00951 TEST_FOR_EXCEPTION( ( ~isValidCompEngine(compEngine) ), std::invalid_argument, 00952 ">>> ERROR (ArrayTools::contractDataDataScalar): Computational engine not defined!"); 00953 } // switch(compEngine) 00954 } // end contractDataDataScalar 00955 00956 00957 template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight> 00958 void ArrayTools::contractDataDataVector(ArrayOutData & outputData, 00959 const ArrayInDataLeft & inputDataLeft, 00960 const ArrayInDataRight & inputDataRight, 00961 const ECompEngine compEngine, 00962 const bool sumInto) { 00963 00964 #ifdef HAVE_INTREPID_DEBUG 00965 TEST_FOR_EXCEPTION( (inputDataLeft.rank() != 3 ), std::invalid_argument, 00966 ">>> ERROR (ArrayTools::contractDataDataVector): Rank of the left input argument must equal 3!"); 00967 TEST_FOR_EXCEPTION( (inputDataRight.rank() != 3 ), std::invalid_argument, 00968 ">>> ERROR (ArrayTools::contractDataDataVector): Rank of right input argument must equal 3!"); 00969 TEST_FOR_EXCEPTION( (outputData.rank() != 1 ), std::invalid_argument, 00970 ">>> ERROR (ArrayTools::contractDataDataVector): Rank of output argument must equal 1!"); 00971 TEST_FOR_EXCEPTION( (inputDataLeft.dimension(0) != inputDataRight.dimension(0) ), std::invalid_argument, 00972 ">>> ERROR (ArrayTools::contractDataDataVector): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!"); 00973 TEST_FOR_EXCEPTION( (inputDataLeft.dimension(1) != inputDataRight.dimension(1) ), std::invalid_argument, 00974 ">>> ERROR (ArrayTools::contractDataDataVector): First dimensions (numbers of integration points) of the left and right input containers must agree!"); 00975 TEST_FOR_EXCEPTION( (inputDataLeft.dimension(2) != inputDataRight.dimension(2) ), std::invalid_argument, 00976 ">>> ERROR (ArrayTools::contractDataDataVector): Second dimensions (numbers of vector components) of the left and right input containers must agree!"); 00977 TEST_FOR_EXCEPTION( (outputData.dimension(0) != inputDataRight.dimension(0) ), std::invalid_argument, 00978 ">>> ERROR (ArrayTools::contractDataDataVector): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!"); 00979 #endif 00980 00981 // get sizes 00982 int numCells = inputDataLeft.dimension(0); 00983 int numPoints = inputDataLeft.dimension(1); 00984 int dimVec = inputDataLeft.dimension(2); 00985 00986 switch(compEngine) { 00987 case COMP_CPP: { 00988 if (sumInto) { 00989 for (int cl = 0; cl < numCells; cl++) { 00990 Scalar tmpVal(0); 00991 for (int qp = 0; qp < numPoints; qp++) { 00992 for (int iVec = 0; iVec < dimVec; iVec++) { 00993 tmpVal += inputDataLeft(cl, qp, iVec)*inputDataRight(cl, qp, iVec); 00994 } // D-loop 00995 } // P-loop 00996 outputData(cl) += tmpVal; 00997 } // C-loop 00998 } 00999 else { 01000 for (int cl = 0; cl < numCells; cl++) { 01001 Scalar tmpVal(0); 01002 for (int qp = 0; qp < numPoints; qp++) { 01003 for (int iVec = 0; iVec < dimVec; iVec++) { 01004 tmpVal += inputDataLeft(cl, qp, iVec)*inputDataRight(cl, qp, iVec); 01005 } // D-loop 01006 } // P-loop 01007 outputData(cl) = tmpVal; 01008 } // C-loop 01009 } 01010 } 01011 break; 01012 01013 case COMP_BLAS: { 01014 int skip = numPoints*dimVec; // size of the left data chunk per cell 01015 int incr = 1; // increment 01016 if (sumInto) { 01017 for (int cl=0; cl < numCells; cl++) { 01018 Teuchos::BLAS<int, Scalar> myblas; 01019 outputData(cl) += myblas.DOT(skip, &inputDataLeft[cl*skip], incr, &inputDataRight[cl*skip], incr); 01020 } 01021 } 01022 else { 01023 for (int cl=0; cl < numCells; cl++) { 01024 Teuchos::BLAS<int, Scalar> myblas; 01025 outputData(cl) = myblas.DOT(skip, &inputDataLeft[cl*skip], incr, &inputDataRight[cl*skip], incr); 01026 } 01027 } 01028 } 01029 break; 01030 01031 default: 01032 TEST_FOR_EXCEPTION( ( ~isValidCompEngine(compEngine) ), std::invalid_argument, 01033 ">>> ERROR (ArrayTools::contractDataDataVector): Computational engine not defined!"); 01034 } // switch(compEngine) 01035 } // end contractDataDataVector 01036 01037 01038 template<class Scalar, class ArrayOutData, class ArrayInDataLeft, class ArrayInDataRight> 01039 void ArrayTools::contractDataDataTensor(ArrayOutData & outputData, 01040 const ArrayInDataLeft & inputDataLeft, 01041 const ArrayInDataRight & inputDataRight, 01042 const ECompEngine compEngine, 01043 const bool sumInto) { 01044 01045 #ifdef HAVE_INTREPID_DEBUG 01046 TEST_FOR_EXCEPTION( (inputDataLeft.rank() != 4 ), std::invalid_argument, 01047 ">>> ERROR (ArrayTools::contractDataDataTensor): Rank of the left input argument must equal 4"); 01048 TEST_FOR_EXCEPTION( (inputDataRight.rank() != 4 ), std::invalid_argument, 01049 ">>> ERROR (ArrayTools::contractDataDataTensor): Rank of right input argument must equal 4!"); 01050 TEST_FOR_EXCEPTION( (outputData.rank() != 1 ), std::invalid_argument, 01051 ">>> ERROR (ArrayTools::contractDataDataTensor): Rank of output argument must equal 1!"); 01052 TEST_FOR_EXCEPTION( (inputDataLeft.dimension(0) != inputDataRight.dimension(0) ), std::invalid_argument, 01053 ">>> ERROR (ArrayTools::contractDataDataTensor): Zeroth dimensions (number of integration domains) of the left and right input containers must agree!"); 01054 TEST_FOR_EXCEPTION( (inputDataLeft.dimension(1) != inputDataRight.dimension(1) ), std::invalid_argument, 01055 ">>> ERROR (ArrayTools::contractDataDataTensor): First dimensions (numbers of integration points) of the left and right input containers must agree!"); 01056 TEST_FOR_EXCEPTION( (inputDataLeft.dimension(2) != inputDataRight.dimension(2) ), std::invalid_argument, 01057 ">>> ERROR (ArrayTools::contractDataDataTensor): Second dimensions (first tensor dimensions) of the left and right input containers must agree!"); 01058 TEST_FOR_EXCEPTION( (inputDataLeft.dimension(3) != inputDataRight.dimension(3) ), std::invalid_argument, 01059 ">>> ERROR (ArrayTools::contractDataDataTensor): Third dimensions (second tensor dimensions) of the left and right input containers must agree!"); 01060 TEST_FOR_EXCEPTION( (outputData.dimension(0) != inputDataRight.dimension(0) ), std::invalid_argument, 01061 ">>> ERROR (ArrayTools::contractDataDataTensor): Zeroth dimensions (numbers of integration domains) of the input and output containers must agree!"); 01062 #endif 01063 01064 // get sizes 01065 int numCells = inputDataLeft.dimension(0); 01066 int numPoints = inputDataLeft.dimension(1); 01067 int dim1Tensor = inputDataLeft.dimension(2); 01068 int dim2Tensor = inputDataLeft.dimension(3); 01069 01070 switch(compEngine) { 01071 case COMP_CPP: { 01072 if (sumInto) { 01073 for (int cl = 0; cl < numCells; cl++) { 01074 Scalar tmpVal(0); 01075 for (int qp = 0; qp < numPoints; qp++) { 01076 for (int iTens1 = 0; iTens1 < dim1Tensor; iTens1++) { 01077 for (int iTens2 = 0; iTens2 < dim2Tensor; iTens2++) { 01078 tmpVal += inputDataLeft(cl, qp, iTens1, iTens2)*inputDataRight(cl, qp, iTens1, iTens2); 01079 } // D2-loop 01080 } // D1-loop 01081 } // P-loop 01082 outputData(cl) += tmpVal; 01083 } // C-loop 01084 } 01085 else { 01086 for (int cl = 0; cl < numCells; cl++) { 01087 Scalar tmpVal(0); 01088 for (int qp = 0; qp < numPoints; qp++) { 01089 for (int iTens1 = 0; iTens1 < dim1Tensor; iTens1++) { 01090 for (int iTens2 = 0; iTens2 < dim2Tensor; iTens2++) { 01091 tmpVal += inputDataLeft(cl, qp, iTens1, iTens2)*inputDataRight(cl, qp, iTens1, iTens2); 01092 } // D2-loop 01093 } // D1-loop 01094 } // P-loop 01095 outputData(cl) = tmpVal; 01096 } // C-loop 01097 } 01098 } 01099 break; 01100 01101 case COMP_BLAS: { 01102 int skip = numPoints*dim1Tensor*dim2Tensor; // size of the left data chunk per cell 01103 int incr = 1; // increment 01104 if (sumInto) { 01105 for (int cl=0; cl < numCells; cl++) { 01106 Teuchos::BLAS<int, Scalar> myblas; 01107 outputData(cl) += myblas.DOT(skip, &inputDataLeft[cl*skip], incr, &inputDataRight[cl*skip], incr); 01108 } 01109 } 01110 else { 01111 for (int cl=0; cl < numCells; cl++) { 01112 Teuchos::BLAS<int, Scalar> myblas; 01113 outputData(cl) = myblas.DOT(skip, &inputDataLeft[cl*skip], incr, &inputDataRight[cl*skip], incr); 01114 } 01115 } 01116 } 01117 break; 01118 01119 default: 01120 TEST_FOR_EXCEPTION( ( ~isValidCompEngine(compEngine) ), std::invalid_argument, 01121 ">>> ERROR (ArrayTools::contractDataDataTensor): Computational engine not defined!"); 01122 } // switch(compEngine) 01123 } // end contractDataDataTensor 01124 01125 01126 } // end namespace Intrepid
1.7.4