Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/src/Shared/Intrepid_ArrayToolsDefContractions.hpp
Go to the documentation of this file.
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