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