Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/src/Shared/Intrepid_RealSpaceToolsDef.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 
00050 namespace Intrepid {
00051 
00052 
00053 
00054 template<class Scalar>
00055 void RealSpaceTools<Scalar>::absval(Scalar* absArray, const Scalar* inArray, const int size) {
00056   for (int i=0; i<size; i++) {
00057     absArray[i] = std::abs(inArray[i]);
00058   }
00059 }
00060 
00061 
00062 
00063 template<class Scalar>
00064 void RealSpaceTools<Scalar>::absval(Scalar* inoutAbsArray, const int size) {
00065   for (int i=0; i<size; i++) {
00066     inoutAbsArray[i] = std::abs(inoutAbsArray[i]);
00067   }
00068 }
00069 
00070 
00071 
00072 template<class Scalar>
00073 template<class ArrayAbs, class ArrayIn>
00074 void RealSpaceTools<Scalar>::absval(ArrayAbs & absArray, const ArrayIn & inArray) {
00075 #ifdef HAVE_INTREPID_DEBUG
00076     TEUCHOS_TEST_FOR_EXCEPTION( ( inArray.rank() != absArray.rank() ),
00077                         std::invalid_argument,
00078                         ">>> ERROR (RealSpaceTools::absval): Array arguments must have identical ranks!");
00079     for (int i=0; i<inArray.rank(); i++) {
00080       TEUCHOS_TEST_FOR_EXCEPTION( ( inArray.dimension(i) != absArray.dimension(i) ),
00081                           std::invalid_argument,
00082                           ">>> ERROR (RealSpaceTools::absval): Dimensions of array arguments do not agree!");
00083     }
00084 #endif
00085 
00086   for (int i=0; i<inArray.size(); i++) {
00087     absArray[i] = std::abs(inArray[i]);
00088   }
00089 }
00090 
00091 
00092 
00093 template<class Scalar>
00094 template<class ArrayInOut>
00095 void RealSpaceTools<Scalar>::absval(ArrayInOut & inoutAbsArray) {
00096   for (int i=0; i<inoutAbsArray.size(); i++) {
00097     inoutAbsArray[i] = std::abs(inoutAbsArray[i]);
00098   }
00099 }
00100 
00101 
00102 
00103 template<class Scalar>
00104 Scalar RealSpaceTools<Scalar>::vectorNorm(const Scalar* inVec, const int dim, const ENorm normType) {
00105   Scalar temp = (Scalar)0;
00106   switch(normType) {
00107     case NORM_TWO:
00108       for(int i = 0; i < dim; i++){
00109         temp += inVec[i]*inVec[i];
00110       }
00111       temp = std::sqrt(temp);
00112       break;
00113     case NORM_INF:
00114       temp = std::abs(inVec[0]);
00115       for(int i = 1; i < dim; i++){
00116         Scalar absData = std::abs(inVec[i]);
00117         if (temp < absData) temp = absData;
00118       }
00119       break;
00120     case NORM_ONE:
00121       for(int i = 0; i < dim; i++){
00122         temp += std::abs(inVec[i]);
00123       }
00124       break;
00125     default:
00126       TEUCHOS_TEST_FOR_EXCEPTION( ( (normType != NORM_TWO) && (normType != NORM_INF) && (normType != NORM_ONE) ),
00127                           std::invalid_argument,
00128                           ">>> ERROR (RealSpaceTools::vectorNorm): Invalid argument normType.");
00129   }
00130   return temp;
00131 }
00132 
00133 
00134 
00135 template<class Scalar>
00136 template<class ArrayIn>
00137 Scalar RealSpaceTools<Scalar>::vectorNorm(const ArrayIn & inVec, const ENorm normType) {
00138 
00139 #ifdef HAVE_INTREPID_DEBUG
00140     TEUCHOS_TEST_FOR_EXCEPTION( ( inVec.rank() != 1 ),
00141                         std::invalid_argument,
00142                         ">>> ERROR (RealSpaceTools::vectorNorm): Vector argument must have rank 1!");
00143 #endif
00144 
00145   int size = inVec.size();
00146 
00147   Scalar temp = (Scalar)0;
00148   switch(normType) {
00149     case NORM_TWO:
00150       for(int i = 0; i < size; i++){
00151         temp += inVec[i]*inVec[i];
00152       }
00153       temp = std::sqrt(temp);
00154       break;
00155     case NORM_INF:
00156       temp = std::abs(inVec[0]);
00157       for(int i = 1; i < size; i++){
00158         Scalar absData = std::abs(inVec[i]);
00159         if (temp < absData) temp = absData;
00160       }
00161       break;
00162     case NORM_ONE:
00163       for(int i = 0; i < size; i++){
00164         temp += std::abs(inVec[i]);
00165       }
00166       break;
00167     default:
00168       TEUCHOS_TEST_FOR_EXCEPTION( ( (normType != NORM_TWO) && (normType != NORM_INF) && (normType != NORM_ONE) ),
00169                           std::invalid_argument,
00170                           ">>> ERROR (RealSpaceTools::vectorNorm): Invalid argument normType.");
00171   }
00172   return temp;
00173 }
00174 
00175 
00176 
00177 template<class Scalar>
00178 template<class ArrayNorm, class ArrayIn>
00179 void RealSpaceTools<Scalar>::vectorNorm(ArrayNorm & normArray, const ArrayIn & inVecs, const ENorm normType) {
00180 
00181   int arrayRank = inVecs.rank();
00182 
00183 #ifdef HAVE_INTREPID_DEBUG
00184     TEUCHOS_TEST_FOR_EXCEPTION( ( arrayRank != normArray.rank()+1 ),
00185                         std::invalid_argument,
00186                         ">>> ERROR (RealSpaceTools::vectorNorm): Ranks of norm and vector array arguments are incompatible!");
00187     TEUCHOS_TEST_FOR_EXCEPTION( ( (arrayRank < 2) || (arrayRank > 3) ),
00188                         std::invalid_argument,
00189                         ">>> ERROR (RealSpaceTools::vectorNorm): Rank of vector array must be 2 or 3!");
00190     for (int i=0; i<arrayRank-1; i++) {
00191       TEUCHOS_TEST_FOR_EXCEPTION( ( inVecs.dimension(i) != normArray.dimension(i) ),
00192                           std::invalid_argument,
00193                           ">>> ERROR (RealSpaceTools::vectorNorm): Dimensions of norm and vector arguments do not agree!");
00194     }
00195 #endif
00196 
00197   int dim_i0 = 1; // first  index dimension (e.g. cell index)
00198   int dim_i1 = 1; // second index dimension (e.g. point index)
00199   int dim    = inVecs.dimension(arrayRank-1); // spatial dimension
00200 
00201   // determine i0 and i1 dimensions
00202   switch(arrayRank) {
00203     case 3:
00204       dim_i0 = inVecs.dimension(0);
00205       dim_i1 = inVecs.dimension(1);
00206       break;
00207     case 2:
00208       dim_i1 = inVecs.dimension(0);
00209       break;
00210   }
00211 
00212   switch(normType) {
00213     case NORM_TWO: {
00214       int offset_i0, offset, normOffset;
00215       for (int i0=0; i0<dim_i0; i0++) {
00216         offset_i0 = i0*dim_i1;
00217         for (int i1=0; i1<dim_i1; i1++) {
00218           offset      = offset_i0 + i1;
00219           normOffset  = offset;
00220           offset     *= dim;
00221           Scalar temp = (Scalar)0;
00222           for(int i = 0; i < dim; i++){
00223             temp += inVecs[offset+i]*inVecs[offset+i];
00224           }
00225           normArray[normOffset] = std::sqrt(temp);
00226         }
00227       }
00228       break;
00229     } // case NORM_TWO
00230 
00231     case NORM_INF: {
00232       int offset_i0, offset, normOffset;
00233       for (int i0=0; i0<dim_i0; i0++) {
00234         offset_i0 = i0*dim_i1;
00235         for (int i1=0; i1<dim_i1; i1++) {
00236           offset      = offset_i0 + i1;
00237           normOffset  = offset;
00238           offset     *= dim;
00239           Scalar temp = (Scalar)0;
00240           temp = std::abs(inVecs[offset]);
00241           for(int i = 1; i < dim; i++){
00242             Scalar absData = std::abs(inVecs[offset+i]);
00243             if (temp < absData) temp = absData;
00244           }
00245           normArray[normOffset] = temp;
00246         }
00247       }
00248       break;
00249     } // case NORM_INF
00250 
00251     case NORM_ONE: {
00252       int offset_i0, offset, normOffset;
00253       for (int i0=0; i0<dim_i0; i0++) {
00254         offset_i0 = i0*dim_i1;
00255         for (int i1=0; i1<dim_i1; i1++) {
00256           offset      = offset_i0 + i1;
00257           normOffset  = offset;
00258           offset     *= dim;
00259           Scalar temp = (Scalar)0;
00260           for(int i = 0; i < dim; i++){
00261             temp += std::abs(inVecs[offset+i]);
00262           }
00263           normArray[normOffset] = temp;
00264         }
00265       }
00266       break;
00267     } // case NORM_ONE
00268 
00269     default:
00270       TEUCHOS_TEST_FOR_EXCEPTION( ( (normType != NORM_TWO) && (normType != NORM_INF) && (normType != NORM_ONE) ),
00271                           std::invalid_argument,
00272                           ">>> ERROR (RealSpaceTools::vectorNorm): Invalid argument normType.");
00273   }
00274 }
00275 
00276 
00277 
00278 template<class Scalar>
00279 void RealSpaceTools<Scalar>::transpose(Scalar* transposeMat, const Scalar* inMat, const int dim) {
00280   for(int i=0; i < dim; i++){
00281     transposeMat[i*dim+i]=inMat[i*dim+i];    // Set diagonal elements
00282     for(int j=i+1; j < dim; j++){
00283       transposeMat[i*dim+j]=inMat[j*dim+i];  // Set off-diagonal elements
00284       transposeMat[j*dim+i]=inMat[i*dim+j];
00285     }
00286   }
00287 }
00288 
00289 
00290 
00291 template<class Scalar>
00292 template<class ArrayTranspose, class ArrayIn>
00293 void RealSpaceTools<Scalar>::transpose(ArrayTranspose & transposeMats, const ArrayIn & inMats) {
00294   int arrayRank = inMats.rank();
00295 
00296 #ifdef HAVE_INTREPID_DEBUG
00297     TEUCHOS_TEST_FOR_EXCEPTION( ( arrayRank != transposeMats.rank() ),
00298                         std::invalid_argument,
00299                         ">>> ERROR (RealSpaceTools::transpose): Matrix array arguments do not have identical ranks!");
00300     TEUCHOS_TEST_FOR_EXCEPTION( ( (arrayRank < 2) || (arrayRank > 4) ),
00301                         std::invalid_argument,
00302                         ">>> ERROR (RealSpaceTools::transpose): Rank of matrix array must be 2, 3, or 4!");
00303     for (int i=0; i<arrayRank; i++) {
00304       TEUCHOS_TEST_FOR_EXCEPTION( ( inMats.dimension(i) != transposeMats.dimension(i) ),
00305                           std::invalid_argument,
00306                           ">>> ERROR (RealSpaceTools::transpose): Dimensions of matrix arguments do not agree!");
00307     }
00308     TEUCHOS_TEST_FOR_EXCEPTION( ( inMats.dimension(arrayRank-2) != inMats.dimension(arrayRank-1) ),
00309                         std::invalid_argument,
00310                         ">>> ERROR (RealSpaceTools::transpose): Matrices are not square!");
00311 #endif
00312 
00313   int dim_i0 = 1; // first  index dimension (e.g. cell index)
00314   int dim_i1 = 1; // second index dimension (e.g. point index)
00315   int dim    = inMats.dimension(arrayRank-2); // spatial dimension
00316 
00317   // determine i0 and i1 dimensions
00318   switch(arrayRank) {
00319     case 4:
00320       dim_i0 = inMats.dimension(0);
00321       dim_i1 = inMats.dimension(1);
00322       break;
00323     case 3:
00324       dim_i1 = inMats.dimension(0);
00325       break;
00326   }
00327 
00328   int offset_i0, offset;
00329 
00330   for (int i0=0; i0<dim_i0; i0++) {
00331     offset_i0 = i0*dim_i1;
00332     for (int i1=0; i1<dim_i1; i1++) {
00333       offset  = offset_i0 + i1;
00334       offset *= (dim*dim);
00335 
00336       for(int i=0; i < dim; i++){
00337         transposeMats[offset+i*dim+i]=inMats[offset+i*dim+i];    // Set diagonal elements
00338         for(int j=i+1; j < dim; j++){
00339           transposeMats[offset+i*dim+j]=inMats[offset+j*dim+i];  // Set off-diagonal elements
00340           transposeMats[offset+j*dim+i]=inMats[offset+i*dim+j];
00341         }
00342       }
00343 
00344     } // i1
00345   } // i0
00346 
00347 }
00348 
00349 
00350 
00351 template<class Scalar>
00352 void RealSpaceTools<Scalar>::inverse(Scalar* inverseMat, const Scalar* inMat, const int dim) {
00353 
00354   switch(dim) {
00355     case 3: {
00356       int i, j, rowID = 0, colID = 0;
00357       int rowperm[3]={0,1,2};
00358       int colperm[3]={0,1,2}; // Complete pivoting
00359       Scalar emax(0);
00360 
00361       for(i=0; i < 3; ++i){
00362         for(j=0; j < 3; ++j){
00363           if( std::abs( inMat[i*3+j] ) >  emax){
00364             rowID = i;  colID = j; emax = std::abs( inMat[i*3+j] );
00365           }
00366         }
00367       }
00368 #ifdef HAVE_INTREPID_DEBUG
00369       TEUCHOS_TEST_FOR_EXCEPTION( ( emax == (Scalar)0 ),
00370                           std::invalid_argument,
00371                           ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
00372 #endif
00373       if( rowID ){
00374         rowperm[0] = rowID;
00375         rowperm[rowID] = 0;
00376       }
00377       if( colID ){
00378         colperm[0] = colID;
00379         colperm[colID] = 0;
00380       }
00381       Scalar B[3][3], S[2][2], Bi[3][3]; // B=rowperm inMat colperm, S=Schur complement(Boo)
00382       for(i=0; i < 3; ++i){
00383         for(j=0; j < 3; ++j){
00384           B[i][j] = inMat[rowperm[i]*3+colperm[j]];
00385         }
00386       }
00387       B[1][0] /= B[0][0]; B[2][0] /= B[0][0];// B(:,0)/=pivot
00388       for(i=0; i < 2; ++i){
00389         for(j=0; j < 2; ++j){
00390           S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1]; // S = B -z*y'
00391         }
00392       }
00393       Scalar detS = S[0][0]*S[1][1]- S[0][1]*S[1][0], Si[2][2];
00394 #ifdef HAVE_INTREPID_DEBUG
00395       TEUCHOS_TEST_FOR_EXCEPTION( ( detS == (Scalar)0 ),
00396                           std::invalid_argument,
00397                           ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!");
00398 #endif
00399 
00400       Si[0][0] =  S[1][1]/detS;                  Si[0][1] = -S[0][1]/detS;
00401       Si[1][0] = -S[1][0]/detS;                  Si[1][1] =  S[0][0]/detS;
00402 
00403       for(j=0; j<2;j++)
00404         Bi[0][j+1] = -( B[0][1]*Si[0][j] + B[0][2]* Si[1][j])/B[0][0];
00405       for(i=0; i<2;i++)
00406         Bi[i+1][0] = -(Si[i][0]*B[1][0] + Si[i][1]*B[2][0]);
00407 
00408       Bi[0][0] =  ((Scalar)1/B[0][0])-Bi[0][1]*B[1][0]-Bi[0][2]*B[2][0];
00409       Bi[1][1] =  Si[0][0];
00410       Bi[1][2] =  Si[0][1];
00411       Bi[2][1] =  Si[1][0];
00412       Bi[2][2] =  Si[1][1];
00413       for(i=0; i < 3; ++i){
00414         for(j=0; j < 3; ++j){
00415           inverseMat[i*3+j] = Bi[colperm[i]][rowperm[j]]; // set inverse
00416         }
00417       }
00418       break;
00419     } // case 3
00420 
00421     case 2: {
00422 
00423       Scalar determinant    = inMat[0]*inMat[3]-inMat[1]*inMat[2];;
00424 #ifdef HAVE_INTREPID_DEBUG
00425       TEUCHOS_TEST_FOR_EXCEPTION( ( (inMat[0]==(Scalar)0) && (inMat[1]==(Scalar)0) &&
00426                             (inMat[2]==(Scalar)0) && (inMat[3]==(Scalar)0) ),
00427                           std::invalid_argument,
00428                           ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
00429       TEUCHOS_TEST_FOR_EXCEPTION( ( determinant == (Scalar)0 ),
00430                           std::invalid_argument,
00431                           ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!");
00432 #endif
00433       inverseMat[0] =   inMat[3] / determinant;
00434       inverseMat[1] = - inMat[1] / determinant;
00435       //
00436       inverseMat[2] = - inMat[2] / determinant;
00437       inverseMat[3] =   inMat[0] / determinant;
00438       break;
00439     } // case 2
00440 
00441     case 1: {
00442 #ifdef HAVE_INTREPID_DEBUG
00443       TEUCHOS_TEST_FOR_EXCEPTION( ( inMat[0] == (Scalar)0 ),
00444                           std::invalid_argument,
00445                           ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
00446 #endif
00447       inverseMat[0] = (Scalar)1 / inMat[0];
00448       break;
00449     } // case 1
00450 
00451   } // switch (dim)
00452 }
00453 
00454 
00455 
00456 template<class Scalar>
00457 template<class ArrayInverse, class ArrayIn>
00458 void RealSpaceTools<Scalar>::inverse(ArrayInverse & inverseMats, const ArrayIn & inMats) {
00459 
00460   int arrayRank = inMats.rank();
00461 
00462 #ifdef HAVE_INTREPID_DEBUG
00463     TEUCHOS_TEST_FOR_EXCEPTION( ( arrayRank != inverseMats.rank() ),
00464                         std::invalid_argument,
00465                         ">>> ERROR (RealSpaceTools::inverse): Matrix array arguments do not have identical ranks!");
00466     TEUCHOS_TEST_FOR_EXCEPTION( ( (arrayRank < 2) || (arrayRank > 4) ),
00467                         std::invalid_argument,
00468                         ">>> ERROR (RealSpaceTools::inverse): Rank of matrix array must be 2, 3, or 4!");
00469     for (int i=0; i<arrayRank; i++) {
00470       TEUCHOS_TEST_FOR_EXCEPTION( ( inMats.dimension(i) != inverseMats.dimension(i) ),
00471                           std::invalid_argument,
00472                           ">>> ERROR (RealSpaceTools::inverse): Dimensions of matrix arguments do not agree!");
00473     }
00474     TEUCHOS_TEST_FOR_EXCEPTION( ( inMats.dimension(arrayRank-2) != inMats.dimension(arrayRank-1) ),
00475                         std::invalid_argument,
00476                         ">>> ERROR (RealSpaceTools::inverse): Matrices are not square!");
00477     TEUCHOS_TEST_FOR_EXCEPTION( ( (inMats.dimension(arrayRank-2) < 1) || (inMats.dimension(arrayRank-2) > 3) ),
00478                         std::invalid_argument,
00479                         ">>> ERROR (RealSpaceTools::inverse): Spatial dimension must be 1, 2, or 3!");
00480 #endif
00481 
00482   int dim_i0 = 1; // first  index dimension (e.g. cell index)
00483   int dim_i1 = 1; // second index dimension (e.g. point index)
00484   int dim    = inMats.dimension(arrayRank-2); // spatial dimension
00485 
00486   // determine i0 and i1 dimensions
00487   switch(arrayRank) {
00488     case 4:
00489       dim_i0 = inMats.dimension(0);
00490       dim_i1 = inMats.dimension(1);
00491       break;
00492     case 3:
00493       dim_i1 = inMats.dimension(0);
00494       break;
00495   }
00496 
00497   switch(dim) {
00498     case 3: {
00499       int offset_i0, offset;
00500 
00501       for (int i0=0; i0<dim_i0; i0++) {
00502         offset_i0 = i0*dim_i1;
00503         for (int i1=0; i1<dim_i1; i1++) {
00504           offset  = offset_i0 + i1;
00505           offset *= 9;
00506 
00507           int i, j, rowID = 0, colID = 0;
00508           int rowperm[3]={0,1,2};
00509           int colperm[3]={0,1,2}; // Complete pivoting
00510           Scalar emax(0);
00511 
00512           for(i=0; i < 3; ++i){
00513             for(j=0; j < 3; ++j){
00514               if( std::abs( inMats[offset+i*3+j] ) >  emax){
00515                 rowID = i;  colID = j; emax = std::abs( inMats[offset+i*3+j] );
00516               }
00517             }
00518           }
00519 #ifdef HAVE_INTREPID_DEBUG
00520 #ifdef HAVE_INTREPID_DEBUG_INF_CHECK
00521           TEUCHOS_TEST_FOR_EXCEPTION( ( emax == (Scalar)0 ),
00522                               std::invalid_argument,
00523                               ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
00524 #endif
00525 #endif
00526           if( rowID ){
00527             rowperm[0] = rowID;
00528             rowperm[rowID] = 0;
00529           }
00530           if( colID ){
00531             colperm[0] = colID;
00532             colperm[colID] = 0;
00533           }
00534           Scalar B[3][3], S[2][2], Bi[3][3]; // B=rowperm inMat colperm, S=Schur complement(Boo)
00535           for(i=0; i < 3; ++i){
00536             for(j=0; j < 3; ++j){
00537               B[i][j] = inMats[offset+rowperm[i]*3+colperm[j]];
00538             }
00539           }
00540           B[1][0] /= B[0][0]; B[2][0] /= B[0][0];// B(:,0)/=pivot
00541           for(i=0; i < 2; ++i){
00542             for(j=0; j < 2; ++j){
00543               S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1]; // S = B -z*y'
00544             }
00545           }
00546           Scalar detS = S[0][0]*S[1][1]- S[0][1]*S[1][0], Si[2][2];
00547 #ifdef HAVE_INTREPID_DEBUG
00548 #ifdef HAVE_INTREPID_DEBUG_INF_CHECK
00549           TEUCHOS_TEST_FOR_EXCEPTION( ( detS == (Scalar)0 ),
00550                               std::invalid_argument,
00551                               ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!");
00552 #endif
00553 #endif
00554 
00555           Si[0][0] =  S[1][1]/detS;                  Si[0][1] = -S[0][1]/detS;
00556           Si[1][0] = -S[1][0]/detS;                  Si[1][1] =  S[0][0]/detS;
00557 
00558           for(j=0; j<2;j++)
00559             Bi[0][j+1] = -( B[0][1]*Si[0][j] + B[0][2]* Si[1][j])/B[0][0];
00560           for(i=0; i<2;i++)
00561             Bi[i+1][0] = -(Si[i][0]*B[1][0] + Si[i][1]*B[2][0]);
00562 
00563           Bi[0][0] =  ((Scalar)1/B[0][0])-Bi[0][1]*B[1][0]-Bi[0][2]*B[2][0];
00564           Bi[1][1] =  Si[0][0];
00565           Bi[1][2] =  Si[0][1];
00566           Bi[2][1] =  Si[1][0];
00567           Bi[2][2] =  Si[1][1];
00568           for(i=0; i < 3; ++i){
00569             for(j=0; j < 3; ++j){
00570               inverseMats[offset+i*3+j] = Bi[colperm[i]][rowperm[j]]; // set inverse
00571             }
00572           }
00573         } // for i1
00574       } // for i0
00575       break;
00576     } // case 3
00577 
00578     case 2: {
00579       int offset_i0, offset;
00580 
00581       for (int i0=0; i0<dim_i0; i0++) {
00582         offset_i0 = i0*dim_i1;
00583         for (int i1=0; i1<dim_i1; i1++) {
00584           offset  = offset_i0 + i1;;
00585           offset *= 4;
00586 
00587           Scalar determinant    = inMats[offset]*inMats[offset+3]-inMats[offset+1]*inMats[offset+2];;
00588 #ifdef HAVE_INTREPID_DEBUG
00589 #ifdef HAVE_INTREPID_DEBUG_INF_CHECK
00590           TEUCHOS_TEST_FOR_EXCEPTION( ( (inMats[offset]==(Scalar)0)   && (inMats[offset+1]==(Scalar)0) &&
00591                                 (inMats[offset+2]==(Scalar)0) && (inMats[offset+3]==(Scalar)0) ),
00592                               std::invalid_argument,
00593                               ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
00594           TEUCHOS_TEST_FOR_EXCEPTION( ( determinant == (Scalar)0 ),
00595                               std::invalid_argument,
00596                               ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!");
00597 #endif
00598 #endif
00599           inverseMats[offset]   = inMats[offset+3] / determinant;
00600           inverseMats[offset+1] = - inMats[offset+1] / determinant;
00601           //
00602           inverseMats[offset+2] = - inMats[offset+2] / determinant;
00603           inverseMats[offset+3] =   inMats[offset] / determinant;
00604         } // for i1
00605       } // for i0
00606       break;
00607     } // case 2
00608 
00609     case 1: {
00610       int offset_i0, offset;
00611 
00612       for (int i0=0; i0<dim_i0; i0++) {
00613         offset_i0 = i0*dim_i1;
00614         for (int i1=0; i1<dim_i1; i1++) {
00615           offset  = offset_i0 + i1;;
00616 #ifdef HAVE_INTREPID_DEBUG
00617 #ifdef HAVE_INTREPID_DEBUG_INF_CHECK
00618           TEUCHOS_TEST_FOR_EXCEPTION( ( inMats[offset] == (Scalar)0 ),
00619                               std::invalid_argument,
00620                               ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
00621 #endif
00622 #endif
00623           inverseMats[offset] = (Scalar)1 / inMats[offset];
00624         } // for i1
00625       } // for i2
00626       break;
00627     } // case 1
00628 
00629   } // switch (dim)
00630 }
00631 
00632 
00633 
00634 template<class Scalar>
00635 Scalar RealSpaceTools<Scalar>::det(const Scalar* inMat, const int dim) {
00636   Scalar determinant(0);
00637 
00638   switch (dim) {
00639     case 3: {
00640       int i,j,rowID = 0;
00641       int colID = 0;
00642       int rowperm[3]={0,1,2};
00643       int colperm[3]={0,1,2}; // Complete pivoting
00644       Scalar emax(0);
00645 
00646       for(i=0; i < 3; ++i){
00647         for(j=0; j < 3; ++j){
00648           if( std::abs( inMat[i*dim+j] ) >  emax){
00649             rowID = i;  colID = j; emax = std::abs( inMat[i*dim+j] );
00650           }
00651         }
00652       }
00653       if( emax > 0 ){
00654         if( rowID ){
00655           rowperm[0] = rowID;
00656           rowperm[rowID] = 0;
00657         }
00658         if( colID ){
00659           colperm[0] = colID;
00660           colperm[colID] = 0;
00661         }
00662         Scalar B[3][3], S[2][2]; // B=rowperm inMat colperm, S=Schur complement(Boo)
00663         for(i=0; i < 3; ++i){
00664           for(j=0; j < 3; ++j){
00665             B[i][j] = inMat[rowperm[i]*dim+colperm[j]];
00666           }
00667         }
00668         B[1][0] /= B[0][0]; B[2][0] /= B[0][0];// B(:,0)/=pivot
00669         for(i=0; i < 2; ++i){
00670           for(j=0; j < 2; ++j){
00671             S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1]; // S = B -z*y'
00672           }
00673         }
00674         determinant = B[0][0] * (S[0][0] * S[1][1] - S[0][1] * S[1][0]); // det(B)
00675         if( rowID ) determinant = -determinant;
00676         if( colID ) determinant = -determinant;
00677       }
00678       break;
00679     } // case 3
00680 
00681     case 2:
00682       determinant = inMat[0]*inMat[3]-
00683                     inMat[1]*inMat[2];
00684       break;
00685 
00686     case 1:
00687       determinant = inMat[0];
00688       break;
00689 
00690     default:
00691       TEUCHOS_TEST_FOR_EXCEPTION( ( (dim != 1) && (dim != 2) && (dim != 3) ),
00692                           std::invalid_argument,
00693                           ">>> ERROR (Matrix): Invalid matrix dimension.");
00694   } // switch (dim)
00695 
00696   return determinant;
00697 }
00698 
00699 
00700 
00701 template<class Scalar>
00702 template<class ArrayIn>
00703 Scalar RealSpaceTools<Scalar>::det(const ArrayIn & inMat) {
00704 
00705 #ifdef HAVE_INTREPID_DEBUG
00706     TEUCHOS_TEST_FOR_EXCEPTION( (inMat.rank() != 2),
00707                         std::invalid_argument,
00708                         ">>> ERROR (RealSpaceTools::det): Rank of matrix argument must be 2!");
00709     TEUCHOS_TEST_FOR_EXCEPTION( ( inMat.dimension(0) != inMat.dimension(1) ),
00710                         std::invalid_argument,
00711                         ">>> ERROR (RealSpaceTools::det): Matrix is not square!");
00712     TEUCHOS_TEST_FOR_EXCEPTION( ( (inMat.dimension(0) < 1) || (inMat.dimension(0) > 3) ),
00713                         std::invalid_argument,
00714                         ">>> ERROR (RealSpaceTools::det): Spatial dimension must be 1, 2, or 3!");
00715 #endif
00716 
00717   int dim = inMat.dimension(0);
00718   Scalar determinant(0);
00719 
00720   switch (dim) {
00721     case 3: {
00722       int i,j,rowID = 0;
00723       int colID = 0;
00724       int rowperm[3]={0,1,2};
00725       int colperm[3]={0,1,2}; // Complete pivoting
00726       Scalar emax(0);
00727 
00728       for(i=0; i < 3; ++i){
00729         for(j=0; j < 3; ++j){
00730           if( std::abs( inMat[i*dim+j] ) >  emax){
00731             rowID = i;  colID = j; emax = std::abs( inMat[i*dim+j] );
00732           }
00733         }
00734       }
00735       if( emax > 0 ){
00736         if( rowID ){
00737           rowperm[0] = rowID;
00738           rowperm[rowID] = 0;
00739         }
00740         if( colID ){
00741           colperm[0] = colID;
00742           colperm[colID] = 0;
00743         }
00744         Scalar B[3][3], S[2][2]; // B=rowperm inMat colperm, S=Schur complement(Boo)
00745         for(i=0; i < 3; ++i){
00746           for(j=0; j < 3; ++j){
00747             B[i][j] = inMat[rowperm[i]*dim+colperm[j]];
00748           }
00749         }
00750         B[1][0] /= B[0][0]; B[2][0] /= B[0][0];// B(:,0)/=pivot
00751         for(i=0; i < 2; ++i){
00752           for(j=0; j < 2; ++j){
00753             S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1]; // S = B -z*y'
00754           }
00755         }
00756         determinant = B[0][0] * (S[0][0] * S[1][1] - S[0][1] * S[1][0]); // det(B)
00757         if( rowID ) determinant = -determinant;
00758         if( colID ) determinant = -determinant;
00759       }
00760       break;
00761     } // case 3
00762 
00763     case 2:
00764       determinant = inMat[0]*inMat[3]-
00765                     inMat[1]*inMat[2];
00766       break;
00767 
00768     case 1:
00769       determinant = inMat[0];
00770       break;
00771 
00772     default:
00773       TEUCHOS_TEST_FOR_EXCEPTION( ( (dim != 1) && (dim != 2) && (dim != 3) ),
00774                           std::invalid_argument,
00775                           ">>> ERROR (Matrix): Invalid matrix dimension.");
00776   } // switch (dim)
00777 
00778   return determinant;
00779 }
00780 
00781 
00782 
00783 
00784 template<class Scalar>
00785 template<class ArrayDet, class ArrayIn>
00786 void RealSpaceTools<Scalar>::det(ArrayDet & detArray, const ArrayIn & inMats) {
00787 
00788   int matArrayRank = inMats.rank();
00789 
00790 #ifdef HAVE_INTREPID_DEBUG
00791     TEUCHOS_TEST_FOR_EXCEPTION( ( matArrayRank != detArray.rank()+2 ),
00792                         std::invalid_argument,
00793                         ">>> ERROR (RealSpaceTools::det): Determinant and matrix array arguments do not have compatible ranks!");
00794     TEUCHOS_TEST_FOR_EXCEPTION( ( (matArrayRank < 3) || (matArrayRank > 4) ),
00795                         std::invalid_argument,
00796                         ">>> ERROR (RealSpaceTools::det): Rank of matrix array must be 3 or 4!");
00797     for (int i=0; i<matArrayRank-2; i++) {
00798       TEUCHOS_TEST_FOR_EXCEPTION( ( inMats.dimension(i) != detArray.dimension(i) ),
00799                           std::invalid_argument,
00800                           ">>> ERROR (RealSpaceTools::det): Dimensions of determinant and matrix array arguments do not agree!");
00801     }
00802     TEUCHOS_TEST_FOR_EXCEPTION( ( inMats.dimension(matArrayRank-2) != inMats.dimension(matArrayRank-1) ),
00803                         std::invalid_argument,
00804                         ">>> ERROR (RealSpaceTools::det): Matrices are not square!");
00805     TEUCHOS_TEST_FOR_EXCEPTION( ( (inMats.dimension(matArrayRank-2) < 1) || (inMats.dimension(matArrayRank-2) > 3) ),
00806                         std::invalid_argument,
00807                         ">>> ERROR (RealSpaceTools::det): Spatial dimension must be 1, 2, or 3!");
00808 #endif
00809 
00810   int dim_i0 = 1; // first  index dimension (e.g. cell index)
00811   int dim_i1 = 1; // second index dimension (e.g. point index)
00812   int dim    = inMats.dimension(matArrayRank-2); // spatial dimension
00813 
00814   // determine i0 and i1 dimensions
00815   switch(matArrayRank) {
00816     case 4:
00817       dim_i0 = inMats.dimension(0);
00818       dim_i1 = inMats.dimension(1);
00819       break;
00820     case 3:
00821       dim_i1 = inMats.dimension(0);
00822       break;
00823   }
00824 
00825   switch(dim) {
00826     case 3: {
00827       int offset_i0, offset, detOffset;
00828 
00829       for (int i0=0; i0<dim_i0; i0++) {
00830         offset_i0 = i0*dim_i1;
00831         for (int i1=0; i1<dim_i1; i1++) {
00832           offset     = offset_i0 + i1;
00833           detOffset  = offset;
00834           offset    *= 9;
00835 
00836           int i,j,rowID = 0;
00837           int colID = 0;
00838           int rowperm[3]={0,1,2};
00839           int colperm[3]={0,1,2}; // Complete pivoting
00840           Scalar emax(0), determinant(0);
00841 
00842           for(i=0; i < 3; ++i){
00843             for(j=0; j < 3; ++j){
00844               if( std::abs( inMats[offset+i*3+j] ) >  emax){
00845                 rowID = i;  colID = j; emax = std::abs( inMats[offset+i*3+j] );
00846               }
00847             }
00848           }
00849           if( emax > 0 ){
00850             if( rowID ){
00851               rowperm[0] = rowID;
00852               rowperm[rowID] = 0;
00853             }
00854             if( colID ){
00855               colperm[0] = colID;
00856               colperm[colID] = 0;
00857             }
00858             Scalar B[3][3], S[2][2]; // B=rowperm inMat colperm, S=Schur complement(Boo)
00859             for(i=0; i < 3; ++i){
00860               for(j=0; j < 3; ++j){
00861                 B[i][j] = inMats[offset+rowperm[i]*3+colperm[j]];
00862               }
00863             }
00864             B[1][0] /= B[0][0]; B[2][0] /= B[0][0];// B(:,0)/=pivot
00865             for(i=0; i < 2; ++i){
00866               for(j=0; j < 2; ++j){
00867                 S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1]; // S = B -z*y'
00868               }
00869             }
00870             determinant = B[0][0] * (S[0][0] * S[1][1] - S[0][1] * S[1][0]); // det(B)
00871             if( rowID ) determinant = -determinant;
00872             if( colID ) determinant = -determinant;
00873           }
00874           detArray[detOffset] = determinant;
00875         } // for i1
00876       } // for i0
00877       break;
00878     } // case 3
00879 
00880     case 2: {
00881       int offset_i0, offset, detOffset;
00882 
00883       for (int i0=0; i0<dim_i0; i0++) {
00884         offset_i0 = i0*dim_i1;
00885         for (int i1=0; i1<dim_i1; i1++) {
00886           offset     = offset_i0 + i1;
00887           detOffset  = offset;
00888           offset    *= 4;
00889 
00890           detArray[detOffset] = inMats[offset]*inMats[offset+3]-inMats[offset+1]*inMats[offset+2];;
00891         } // for i1
00892       } // for i0
00893       break;
00894     } // case 2
00895 
00896     case 1: {
00897       int offset_i0, offset;
00898 
00899       for (int i0=0; i0<dim_i0; i0++) {
00900         offset_i0 = i0*dim_i1;
00901         for (int i1=0; i1<dim_i1; i1++) {
00902           offset  = offset_i0 + i1;;
00903           detArray[offset] = inMats[offset];
00904         } // for i1
00905       } // for i2
00906       break;
00907     } // case 1
00908 
00909   } // switch (dim)
00910 }
00911 
00912 
00913 
00914 template<class Scalar>
00915 void RealSpaceTools<Scalar>::add(Scalar* sumArray, const Scalar* inArray1, const Scalar* inArray2, const int size) {
00916   for (int i=0; i<size; i++) {
00917     sumArray[i] = inArray1[i] + inArray2[i];
00918   }
00919 }
00920 
00921 
00922 
00923 template<class Scalar>
00924 void RealSpaceTools<Scalar>::add(Scalar* inoutSumArray, const Scalar* inArray, const int size) {
00925   for (int i=0; i<size; i++) {
00926     inoutSumArray[i] += inArray[i];
00927   }
00928 }
00929 
00930 
00931 
00932 template<class Scalar>
00933 template<class ArraySum, class ArrayIn1, class ArrayIn2>
00934 void RealSpaceTools<Scalar>::add(ArraySum & sumArray, const ArrayIn1 & inArray1, const ArrayIn2 & inArray2) {
00935 #ifdef HAVE_INTREPID_DEBUG
00936     TEUCHOS_TEST_FOR_EXCEPTION( ( (inArray1.rank() != inArray2.rank()) || (inArray1.rank() != sumArray.rank()) ),
00937                         std::invalid_argument,
00938                         ">>> ERROR (RealSpaceTools::add): Array arguments must have identical ranks!");
00939     for (int i=0; i<inArray1.rank(); i++) {
00940       TEUCHOS_TEST_FOR_EXCEPTION( ( (inArray1.dimension(i) != inArray2.dimension(i)) || (inArray1.dimension(i) != sumArray.dimension(i)) ),
00941                           std::invalid_argument,
00942                           ">>> ERROR (RealSpaceTools::add): Dimensions of array arguments do not agree!");
00943     }
00944 #endif
00945 
00946   for (int i=0; i<inArray1.size(); i++) {
00947     sumArray[i] = inArray1[i] + inArray2[i];
00948   }
00949 }
00950 
00951 
00952 
00953 template<class Scalar>
00954 template<class ArraySum, class ArrayIn>
00955 void RealSpaceTools<Scalar>::add(ArraySum & inoutSumArray, const ArrayIn & inArray) {
00956 #ifdef HAVE_INTREPID_DEBUG
00957     TEUCHOS_TEST_FOR_EXCEPTION( ( inArray.rank() != inoutSumArray.rank() ),
00958                         std::invalid_argument,
00959                         ">>> ERROR (RealSpaceTools::add): Array arguments must have identical ranks!");
00960     for (int i=0; i<inArray.rank(); i++) {
00961       TEUCHOS_TEST_FOR_EXCEPTION( ( inArray.dimension(i) != inoutSumArray.dimension(i) ),
00962                           std::invalid_argument,
00963                           ">>> ERROR (RealSpaceTools::add): Dimensions of array arguments do not agree!");
00964     }
00965 #endif
00966 
00967   for (int i=0; i<inArray.size(); i++) {
00968     inoutSumArray[i] += inArray[i];
00969   }
00970 }
00971 
00972 
00973 
00974 template<class Scalar>
00975 void RealSpaceTools<Scalar>::subtract(Scalar* diffArray, const Scalar* inArray1, const Scalar* inArray2, const int size) {
00976   for (int i=0; i<size; i++) {
00977     diffArray[i] = inArray1[i] - inArray2[i];
00978   }
00979 }
00980 
00981 
00982 
00983 template<class Scalar>
00984 void RealSpaceTools<Scalar>::subtract(Scalar* inoutDiffArray, const Scalar* inArray, const int size) {
00985   for (int i=0; i<size; i++) {
00986     inoutDiffArray[i] -= inArray[i];
00987   }
00988 }
00989 
00990 
00991 
00992 template<class Scalar>
00993 template<class ArrayDiff, class ArrayIn1, class ArrayIn2>
00994 void RealSpaceTools<Scalar>::subtract(ArrayDiff & diffArray, const ArrayIn1 & inArray1, const ArrayIn2 & inArray2) {
00995 #ifdef HAVE_INTREPID_DEBUG
00996     TEUCHOS_TEST_FOR_EXCEPTION( ( (inArray1.rank() != inArray2.rank()) || (inArray1.rank() != diffArray.rank()) ),
00997                         std::invalid_argument,
00998                         ">>> ERROR (RealSpaceTools::subtract): Array arguments must have identical ranks!");
00999     for (int i=0; i<inArray1.rank(); i++) {
01000       TEUCHOS_TEST_FOR_EXCEPTION( ( (inArray1.dimension(i) != inArray2.dimension(i)) || (inArray1.dimension(i) != diffArray.dimension(i)) ),
01001                           std::invalid_argument,
01002                           ">>> ERROR (RealSpaceTools::subtract): Dimensions of array arguments do not agree!");
01003     }
01004 #endif
01005 
01006   for (int i=0; i<inArray1.size(); i++) {
01007     diffArray[i] = inArray1[i] - inArray2[i];
01008   }
01009 }
01010 
01011 
01012 
01013 template<class Scalar>
01014 template<class ArrayDiff, class ArrayIn>
01015 void RealSpaceTools<Scalar>::subtract(ArrayDiff & inoutDiffArray, const ArrayIn & inArray) {
01016 #ifdef HAVE_INTREPID_DEBUG
01017     TEUCHOS_TEST_FOR_EXCEPTION( ( inArray.rank() != inoutDiffArray.rank() ),
01018                         std::invalid_argument,
01019                         ">>> ERROR (RealSpaceTools::subtract): Array arguments must have identical ranks!");
01020     for (int i=0; i<inArray.rank(); i++) {
01021       TEUCHOS_TEST_FOR_EXCEPTION( ( inArray.dimension(i) != inoutDiffArray.dimension(i) ),
01022                           std::invalid_argument,
01023                           ">>> ERROR (RealSpaceTools::subtract): Dimensions of array arguments do not agree!");
01024     }
01025 #endif
01026 
01027   for (int i=0; i<inArray.size(); i++) {
01028     inoutDiffArray[i] -= inArray[i];
01029   }
01030 }
01031 
01032 
01033 
01034 
01035 template<class Scalar>
01036 void RealSpaceTools<Scalar>::scale(Scalar* scaledArray, const Scalar* inArray, const int size, const Scalar scalar) {
01037   for (int i=0; i<size; i++) {
01038     scaledArray[i] = scalar*inArray[i];
01039   }
01040 }
01041 
01042 
01043 
01044 template<class Scalar>
01045 void RealSpaceTools<Scalar>::scale(Scalar* inoutScaledArray, const int size, const Scalar scalar) {
01046   for (int i=0; i<size; i++) {
01047     inoutScaledArray[i] *= scalar;
01048   }
01049 }
01050 
01051 
01052 
01053 template<class Scalar>
01054 template<class ArrayScaled, class ArrayIn>
01055 void RealSpaceTools<Scalar>::scale(ArrayScaled & scaledArray, const ArrayIn & inArray, const Scalar scalar) {
01056 #ifdef HAVE_INTREPID_DEBUG
01057     TEUCHOS_TEST_FOR_EXCEPTION( ( inArray.rank() != scaledArray.rank() ),
01058                         std::invalid_argument,
01059                         ">>> ERROR (RealSpaceTools::scale): Array arguments must have identical ranks!");
01060     for (int i=0; i<inArray.rank(); i++) {
01061       TEUCHOS_TEST_FOR_EXCEPTION( ( inArray.dimension(i) != scaledArray.dimension(i) ),
01062                           std::invalid_argument,
01063                           ">>> ERROR (RealSpaceTools::scale): Dimensions of array arguments do not agree!");
01064     }
01065 #endif
01066 
01067   for (int i=0; i<inArray.size(); i++) {
01068     scaledArray[i] = scalar*inArray[i];
01069   }
01070 }
01071 
01072 
01073 
01074 template<class Scalar>
01075 template<class ArrayScaled>
01076 void RealSpaceTools<Scalar>::scale(ArrayScaled & inoutScaledArray, const Scalar scalar) {
01077   for (int i=0; i<inoutScaledArray.size(); i++) {
01078     inoutScaledArray[i] *= scalar;
01079   }
01080 }
01081 
01082 
01083 
01084 
01085 template<class Scalar>
01086 Scalar RealSpaceTools<Scalar>::dot(const Scalar* inArray1, const Scalar* inArray2, const int size) {
01087   Scalar dot(0);
01088   for (int i=0; i<size; i++) {
01089     dot += inArray1[i]*inArray2[i];
01090   }
01091   return dot;  
01092 }
01093 
01094 
01095 
01096 template<class Scalar>
01097 template<class ArrayVec1, class ArrayVec2>
01098 Scalar RealSpaceTools<Scalar>::dot(const ArrayVec1 & inVec1, const ArrayVec2 & inVec2) {
01099 #ifdef HAVE_INTREPID_DEBUG
01100     TEUCHOS_TEST_FOR_EXCEPTION( ( (inVec1.rank() != 1) || (inVec2.rank() != 1) ),
01101                         std::invalid_argument,
01102                         ">>> ERROR (RealSpaceTools::dot): Vector arguments must have rank 1!");
01103     TEUCHOS_TEST_FOR_EXCEPTION( ( inVec1.dimension(0) != inVec2.dimension(0) ),
01104                         std::invalid_argument,
01105                         ">>> ERROR (RealSpaceTools::dot): Dimensions of vector arguments must agree!");
01106 #endif
01107 
01108   Scalar dot(0);
01109   for (int i=0; i<inVec1.size(); i++) {
01110     dot += inVec1[i]*inVec2[i];
01111   }
01112   return dot;  
01113 
01114 }
01115 
01116 
01117 
01118 template<class Scalar>
01119 template<class ArrayDot, class ArrayVec1, class ArrayVec2>
01120 void RealSpaceTools<Scalar>::dot(ArrayDot & dotArray, const ArrayVec1 & inVecs1, const ArrayVec2 & inVecs2) {
01121 
01122   int arrayRank = inVecs1.rank();
01123 
01124 #ifdef HAVE_INTREPID_DEBUG
01125     TEUCHOS_TEST_FOR_EXCEPTION( ( arrayRank != dotArray.rank()+1 ),
01126                         std::invalid_argument,
01127                         ">>> ERROR (RealSpaceTools::dot): Ranks of norm and vector array arguments are incompatible!");
01128     TEUCHOS_TEST_FOR_EXCEPTION( ( arrayRank != inVecs2.rank() ),
01129                         std::invalid_argument,
01130                         ">>> ERROR (RealSpaceTools::dot): Ranks of input vector arguments must be identical!");
01131     TEUCHOS_TEST_FOR_EXCEPTION( ( (arrayRank < 2) || (arrayRank > 3) ),
01132                         std::invalid_argument,
01133                         ">>> ERROR (RealSpaceTools::dot): Rank of input vector arguments must be 2 or 3!");
01134     for (int i=0; i<arrayRank; i++) {
01135       TEUCHOS_TEST_FOR_EXCEPTION( ( inVecs1.dimension(i) != inVecs2.dimension(i) ),
01136                           std::invalid_argument,
01137                           ">>> ERROR (RealSpaceTools::dot): Dimensions of input vector arguments do not agree!");
01138     }
01139     for (int i=0; i<arrayRank-1; i++) {
01140       TEUCHOS_TEST_FOR_EXCEPTION( ( inVecs1.dimension(i) != dotArray.dimension(i) ),
01141                           std::invalid_argument,
01142                           ">>> ERROR (RealSpaceTools::dot): Dimensions of dot-product and vector arrays do not agree!");
01143     }
01144 #endif
01145 
01146   int dim_i0 = 1; // first  index dimension (e.g. cell index)
01147   int dim_i1 = 1; // second index dimension (e.g. point index)
01148   int dim    = inVecs1.dimension(arrayRank-1); // spatial dimension
01149 
01150   // determine i0 and i1 dimensions
01151   switch(arrayRank) {
01152     case 3:
01153       dim_i0 = inVecs1.dimension(0);
01154       dim_i1 = inVecs1.dimension(1);
01155       break;
01156     case 2:
01157       dim_i1 = inVecs1.dimension(0);
01158       break;
01159   }
01160 
01161   int offset_i0, offset, dotOffset;
01162   for (int i0=0; i0<dim_i0; i0++) {
01163     offset_i0 = i0*dim_i1;
01164     for (int i1=0; i1<dim_i1; i1++) {
01165       offset      = offset_i0 + i1;
01166       dotOffset   = offset;
01167       offset     *= dim;
01168       Scalar dot(0);
01169       for (int i=0; i<dim; i++) {
01170         dot += inVecs1[offset+i]*inVecs2[offset+i];
01171       }
01172       dotArray[dotOffset] = dot;
01173     }
01174   }
01175 }
01176 
01177 
01178 
01179 template<class Scalar>
01180 void RealSpaceTools<Scalar>::matvec(Scalar* matVec, const Scalar* inMat, const Scalar* inVec, const int dim) {
01181   for (int i=0; i<dim; i++) {
01182     Scalar sumdot(0);
01183     for (int j=0; j<dim; j++) {
01184       sumdot += inMat[i*dim+j]*inVec[j];
01185     }
01186     matVec[i] = sumdot; 
01187   }
01188 }
01189 
01190 
01191 
01192 template<class Scalar>
01193 template<class ArrayMatVec, class ArrayMat, class ArrayVec>
01194 void RealSpaceTools<Scalar>::matvec(ArrayMatVec & matVecs, const ArrayMat & inMats, const ArrayVec & inVecs) {
01195   int matArrayRank = inMats.rank();
01196 
01197 #ifdef HAVE_INTREPID_DEBUG
01198     TEUCHOS_TEST_FOR_EXCEPTION( ( matArrayRank != inVecs.rank()+1 ),
01199                         std::invalid_argument,
01200                         ">>> ERROR (RealSpaceTools::matvec): Vector and matrix array arguments do not have compatible ranks!");
01201     TEUCHOS_TEST_FOR_EXCEPTION( ( (matArrayRank < 3) || (matArrayRank > 4) ),
01202                         std::invalid_argument,
01203                         ">>> ERROR (RealSpaceTools::matvec): Rank of matrix array must be 3 or 4!");
01204     TEUCHOS_TEST_FOR_EXCEPTION( ( matVecs.rank() != inVecs.rank() ),
01205                         std::invalid_argument,
01206                         ">>> ERROR (RealSpaceTools::matvec): Vector arrays must be have the same rank!");
01207     for (int i=0; i<matArrayRank-1; i++) {
01208       TEUCHOS_TEST_FOR_EXCEPTION( ( inMats.dimension(i) != inVecs.dimension(i) ),
01209                           std::invalid_argument,
01210                           ">>> ERROR (RealSpaceTools::matvec): Dimensions of vector and matrix array arguments do not agree!");
01211     }
01212     for (int i=0; i<inVecs.rank(); i++) {
01213       TEUCHOS_TEST_FOR_EXCEPTION( ( matVecs.dimension(i) != inVecs.dimension(i) ),
01214                           std::invalid_argument,
01215                           ">>> ERROR (RealSpaceTools::matvec): Dimensions of vector array arguments do not agree!");
01216     }
01217     TEUCHOS_TEST_FOR_EXCEPTION( ( inMats.dimension(matArrayRank-2) != inMats.dimension(matArrayRank-1) ),
01218                         std::invalid_argument,
01219                         ">>> ERROR (RealSpaceTools::matvec): Matrices are not square!");
01220 #endif
01221 
01222   int dim_i0 = 1; // first  index dimension (e.g. cell index)
01223   int dim_i1 = 1; // second index dimension (e.g. point index)
01224   int dim    = inMats.dimension(matArrayRank-2); // spatial dimension
01225 
01226   // determine i0 and i1 dimensions
01227   switch(matArrayRank) {
01228     case 4:
01229       dim_i0 = inMats.dimension(0);
01230       dim_i1 = inMats.dimension(1);
01231       break;
01232     case 3:
01233       dim_i1 = inMats.dimension(0);
01234       break;
01235   }
01236 
01237   int offset_i0, offset, vecOffset;
01238 
01239   for (int i0=0; i0<dim_i0; i0++) {
01240     offset_i0 = i0*dim_i1;
01241     for (int i1=0; i1<dim_i1; i1++) {
01242       offset     = offset_i0 + i1;
01243       vecOffset  = offset*dim;
01244       offset     = vecOffset*dim;
01245 
01246       for (int i=0; i<dim; i++) {
01247         Scalar sumdot(0);
01248         for (int j=0; j<dim; j++) {
01249           sumdot += inMats[offset+i*dim+j]*inVecs[vecOffset+j];
01250         }
01251         matVecs[vecOffset+i] = sumdot;
01252       }
01253     }
01254   }
01255 }
01256 
01257 
01258 template<class Scalar>
01259 template<class ArrayVecProd, class ArrayIn1, class ArrayIn2>
01260 void RealSpaceTools<Scalar>::vecprod(ArrayVecProd & vecProd, const ArrayIn1 & inLeft, const ArrayIn2 & inRight) {
01261   
01262 #ifdef HAVE_INTREPID_DEBUG
01263   /*
01264    *   Check array rank and spatial dimension range (if applicable)
01265    *      (1) all array arguments are required to have matching dimensions and rank: (D), (I0,D) or (I0,I1,D)
01266    *      (2) spatial dimension should be 2 or 3
01267    */
01268   std::string errmsg = ">>> ERROR (RealSpaceTools::vecprod):";
01269   
01270   // (1) check rank range on inLeft and then compare the other arrays with inLeft
01271   TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inLeft,  1,3), std::invalid_argument, errmsg);
01272   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inLeft, inRight), std::invalid_argument, errmsg);    
01273   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inLeft, vecProd), std::invalid_argument, errmsg);   
01274   
01275   // (2) spatial dimension ordinal = array rank - 1. Suffices to check only one array because we just
01276   //     checked whether or not the arrays have matching dimensions. 
01277   TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inLeft, inLeft.rank() - 1,  2,3), std::invalid_argument, errmsg);
01278   
01279 #endif
01280 
01281  int spaceDim = inLeft.dimension(inLeft.rank() - 1);
01282 
01283   switch(inLeft.rank() ){
01284     
01285     case 1:
01286       {        
01287         vecProd(0) = inLeft(1)*inRight(2) - inLeft(2)*inRight(1);
01288         vecProd(1) = inLeft(2)*inRight(0) - inLeft(0)*inRight(2);              
01289         vecProd(2) = inLeft(0)*inRight(1) - inLeft(1)*inRight(0);    
01290       }
01291       break;
01292       
01293     case 2:
01294       {
01295         int dim0 = inLeft.dimension(0);
01296         if(spaceDim == 3) {
01297           for(int i0 = 0; i0 < dim0; i0++){
01298             vecProd(i0, 0) = inLeft(i0, 1)*inRight(i0, 2) - inLeft(i0, 2)*inRight(i0, 1);
01299             vecProd(i0, 1) = inLeft(i0, 2)*inRight(i0, 0) - inLeft(i0, 0)*inRight(i0, 2);              
01300             vecProd(i0, 2) = inLeft(i0, 0)*inRight(i0, 1) - inLeft(i0, 1)*inRight(i0, 0);
01301           }// i0
01302         } //spaceDim == 3
01303         else if(spaceDim == 2){
01304           for(int i0 = 0; i0 < dim0; i0++){
01305             // vecprod is scalar - do we still want result to be (i0,i1,D)?
01306             vecProd(i0, 0) = inLeft(i0, 0)*inRight(i0, 1) - inLeft(i0, 1)*inRight(i0, 0);
01307           }// i0
01308         }// spaceDim == 2
01309       }// case 2
01310       break;
01311       
01312     case 3:
01313       {
01314         int dim0 = inLeft.dimension(0);
01315         int dim1 = inLeft.dimension(1);
01316         if(spaceDim == 3) {
01317           for(int i0 = 0; i0 < dim0; i0++){
01318             for(int i1 = 0; i1 < dim1; i1++){
01319               vecProd(i0, i1, 0) = inLeft(i0, i1, 1)*inRight(i0, i1, 2) - inLeft(i0, i1, 2)*inRight(i0, i1, 1);
01320               vecProd(i0, i1, 1) = inLeft(i0, i1, 2)*inRight(i0, i1, 0) - inLeft(i0, i1, 0)*inRight(i0, i1, 2);              
01321               vecProd(i0, i1, 2) = inLeft(i0, i1, 0)*inRight(i0, i1, 1) - inLeft(i0, i1, 1)*inRight(i0, i1, 0);
01322             }// i1
01323           }// i0
01324         } //spaceDim == 3
01325         else if(spaceDim == 2){
01326           for(int i0 = 0; i0 < dim0; i0++){
01327             for(int i1 = 0; i1 < dim1; i1++){
01328               // vecprod is scalar - do we still want result to be (i0,i1,D)?
01329               vecProd(i0, i1, 0) = inLeft(i0, i1, 0)*inRight(i0, i1, 1) - inLeft(i0, i1, 1)*inRight(i0, i1, 0);
01330             }// i1
01331           }// i0
01332         }// spaceDim == 2
01333       } // case 3
01334       break;
01335       
01336     default:
01337       TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, 
01338                          ">>> ERROR (RealSpaceTools::vecprod): rank-1,2,3 arrays required");      
01339   }
01340   
01341 }
01342 
01343 } // namespace Intrepid