00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
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 ArrayScalar>
00060 void RealSpaceTools<Scalar>::absval(ArrayScalar & absArray, const ArrayScalar & 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 ArrayScalar>
00081 void RealSpaceTools<Scalar>::absval(ArrayScalar & 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 VecArray>
00123 Scalar RealSpaceTools<Scalar>::vectorNorm(const VecArray & 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 NormArray, class VecArray>
00165 void RealSpaceTools<Scalar>::vectorNorm(NormArray & normArray, const VecArray & 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;
00184 int dim_i1 = 1;
00185 int dim = inVecs.dimension(arrayRank-1);
00186
00187
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 }
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 }
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 }
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];
00268 for(int j=i+1; j < dim; j++){
00269 transposeMat[i*dim+j]=inMat[j*dim+i];
00270 transposeMat[j*dim+i]=inMat[i*dim+j];
00271 }
00272 }
00273 }
00274
00275
00276
00277 template<class Scalar>
00278 template<class MatArray>
00279 void RealSpaceTools<Scalar>::transpose(MatArray & transposeMats, const MatArray & 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;
00300 int dim_i1 = 1;
00301 int dim = inMats.dimension(arrayRank-2);
00302
00303
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];
00324 for(int j=i+1; j < dim; j++){
00325 transposeMats[offset+i*dim+j]=inMats[offset+j*dim+i];
00326 transposeMats[offset+j*dim+i]=inMats[offset+i*dim+j];
00327 }
00328 }
00329
00330 }
00331 }
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};
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];
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];
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];
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]];
00402 }
00403 }
00404 break;
00405 }
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 }
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 }
00436
00437 }
00438 }
00439
00440
00441
00442 template<class Scalar>
00443 template<class MatArray>
00444 void RealSpaceTools<Scalar>::inverse(MatArray & inverseMats, const MatArray & 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;
00469 int dim_i1 = 1;
00470 int dim = inMats.dimension(arrayRank-2);
00471
00472
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};
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];
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];
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];
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]];
00553 }
00554 }
00555 }
00556 }
00557 break;
00558 }
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 }
00585 }
00586 break;
00587 }
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 }
00603 }
00604 break;
00605 }
00606
00607 }
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};
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];
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];
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];
00650 }
00651 }
00652 determinant = B[0][0] * (S[0][0] * S[1][1] - S[0][1] * S[1][0]);
00653 if( rowID ) determinant = -determinant;
00654 if( colID ) determinant = -determinant;
00655 }
00656 break;
00657 }
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 }
00673
00674 return determinant;
00675 }
00676
00677
00678
00679 template<class Scalar>
00680 template<class ArrayScalar>
00681 Scalar RealSpaceTools<Scalar>::det(const ArrayScalar & 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};
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];
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];
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];
00732 }
00733 }
00734 determinant = B[0][0] * (S[0][0] * S[1][1] - S[0][1] * S[1][0]);
00735 if( rowID ) determinant = -determinant;
00736 if( colID ) determinant = -determinant;
00737 }
00738 break;
00739 }
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 }
00755
00756 return determinant;
00757 }
00758
00759
00760
00761
00762 template<class Scalar>
00763 template<class DetArray, class MatArray>
00764 void RealSpaceTools<Scalar>::det(DetArray & detArray, const MatArray & 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;
00789 int dim_i1 = 1;
00790 int dim = inMats.dimension(matArrayRank-2);
00791
00792
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};
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];
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];
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];
00846 }
00847 }
00848 determinant = B[0][0] * (S[0][0] * S[1][1] - S[0][1] * S[1][0]);
00849 if( rowID ) determinant = -determinant;
00850 if( colID ) determinant = -determinant;
00851 }
00852 detArray[detOffset] = determinant;
00853 }
00854 }
00855 break;
00856 }
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 }
00870 }
00871 break;
00872 }
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 }
00883 }
00884 break;
00885 }
00886
00887 }
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 ArrayScalar>
00912 void RealSpaceTools<Scalar>::add(ArrayScalar & sumArray, const ArrayScalar & inArray1, const ArrayScalar & 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 ArrayScalar>
00933 void RealSpaceTools<Scalar>::add(ArrayScalar & inoutSumArray, const ArrayScalar & 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 ArrayScalar>
00972 void RealSpaceTools<Scalar>::subtract(ArrayScalar & diffArray, const ArrayScalar & inArray1, const ArrayScalar & 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 ArrayScalar>
00993 void RealSpaceTools<Scalar>::subtract(ArrayScalar & inoutDiffArray, const ArrayScalar & 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 ArrayScalar>
01033 void RealSpaceTools<Scalar>::scale(ArrayScalar & scaledArray, const ArrayScalar & 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 ArrayScalar>
01054 void RealSpaceTools<Scalar>::scale(ArrayScalar & 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 VecArray>
01076 Scalar RealSpaceTools<Scalar>::dot(const VecArray & inVec1, const VecArray & 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 DotArray, class VecArray>
01098 void RealSpaceTools<Scalar>::dot(DotArray & dotArray, const VecArray & inVecs1, const VecArray & 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;
01125 int dim_i1 = 1;
01126 int dim = inVecs1.dimension(arrayRank-1);
01127
01128
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 MatArray, class VecArray>
01172 void RealSpaceTools<Scalar>::matvec(VecArray & matVecs, const MatArray & inMats, const VecArray & 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;
01201 int dim_i1 = 1;
01202 int dim = inMats.dimension(matArrayRank-2);
01203
01204
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 VecArrayOut, class VecArrayIn>
01238 void RealSpaceTools<Scalar>::vecprod(VecArrayOut & vecProd,
01239 const VecArrayIn & inLeft,
01240 const VecArrayIn & inRight){
01241
01242 #ifdef HAVE_INTREPID_DEBUG
01243
01244
01245
01246
01247
01248 std::string errmsg = ">>> ERROR (RealSpaceTools::vecprod):";
01249
01250
01251 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inLeft, 1,3), std::invalid_argument, errmsg);
01252 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inLeft, inRight), std::invalid_argument, errmsg);
01253 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inLeft, vecProd), std::invalid_argument, errmsg);
01254
01255
01256
01257 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inLeft, inLeft.rank() - 1, 2,3), std::invalid_argument, errmsg);
01258
01259 #endif
01260
01261 int spaceDim = inLeft.dimension(inLeft.rank() - 1);
01262
01263 switch(inLeft.rank() ){
01264
01265 case 1:
01266 {
01267 vecProd(0) = inLeft(1)*inRight(2) - inLeft(2)*inRight(1);
01268 vecProd(1) = inLeft(2)*inRight(0) - inLeft(0)*inRight(2);
01269 vecProd(2) = inLeft(0)*inRight(1) - inLeft(1)*inRight(0);
01270 }
01271 break;
01272
01273 case 2:
01274 {
01275 int dim0 = inLeft.dimension(0);
01276 if(spaceDim == 3) {
01277 for(int i0 = 0; i0 < dim0; i0++){
01278 vecProd(i0, 0) = inLeft(i0, 1)*inRight(i0, 2) - inLeft(i0, 2)*inRight(i0, 1);
01279 vecProd(i0, 1) = inLeft(i0, 2)*inRight(i0, 0) - inLeft(i0, 0)*inRight(i0, 2);
01280 vecProd(i0, 2) = inLeft(i0, 0)*inRight(i0, 1) - inLeft(i0, 1)*inRight(i0, 0);
01281 }
01282 }
01283 else if(spaceDim == 2){
01284 for(int i0 = 0; i0 < dim0; i0++){
01285
01286 vecProd(i0, 0) = inLeft(i0, 0)*inRight(i0, 1) - inLeft(i0, 1)*inRight(i0, 0);
01287 }
01288 }
01289 }
01290 break;
01291
01292 case 3:
01293 {
01294 int dim0 = inLeft.dimension(0);
01295 int dim1 = inLeft.dimension(1);
01296 if(spaceDim == 3) {
01297 for(int i0 = 0; i0 < dim0; i0++){
01298 for(int i1 = 0; i1 < dim1; i1++){
01299 vecProd(i0, i1, 0) = inLeft(i0, i1, 1)*inRight(i0, i1, 2) - inLeft(i0, i1, 2)*inRight(i0, i1, 1);
01300 vecProd(i0, i1, 1) = inLeft(i0, i1, 2)*inRight(i0, i1, 0) - inLeft(i0, i1, 0)*inRight(i0, i1, 2);
01301 vecProd(i0, i1, 2) = inLeft(i0, i1, 0)*inRight(i0, i1, 1) - inLeft(i0, i1, 1)*inRight(i0, i1, 0);
01302 }
01303 }
01304 }
01305 else if(spaceDim == 2){
01306 for(int i0 = 0; i0 < dim0; i0++){
01307 for(int i1 = 0; i1 < dim1; i1++){
01308
01309 vecProd(i0, i1, 0) = inLeft(i0, i1, 0)*inRight(i0, i1, 1) - inLeft(i0, i1, 1)*inRight(i0, i1, 0);
01310 }
01311 }
01312 }
01313 }
01314 break;
01315
01316 default:
01317 TEST_FOR_EXCEPTION(true, std::invalid_argument,
01318 ">>> ERROR (RealSpaceTools::vecprod): rank-1,2,3 arrays required");
01319 }
01320
01321 }
01322
01323
01324
01325 }
