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