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