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
00030
00031 #include "SundanceProductEvaluator.hpp"
00032 #include "SundanceEvalManager.hpp"
00033 #include "SundanceProductExpr.hpp"
00034
00035 #include "PlayaTabs.hpp"
00036 #include "SundanceOut.hpp"
00037
00038 using namespace Sundance;
00039 using namespace Sundance;
00040 using namespace Sundance;
00041 using namespace Teuchos;
00042
00043
00044
00045
00046
00047 ProductEvaluator::ProductEvaluator(const ProductExpr* expr,
00048 const EvalContext& context)
00049 : BinaryEvaluator<ProductExpr>(expr, context),
00050 maxOrder_(this->sparsity()->maxOrder()),
00051 resultIndex_(maxOrder_+1),
00052 resultIsConstant_(maxOrder_+1),
00053 hasWorkspace_(maxOrder_+1),
00054 workspaceIsLeft_(maxOrder_+1),
00055 workspaceIndex_(maxOrder_+1),
00056 workspaceCoeffIndex_(maxOrder_+1),
00057 workspaceCoeffIsConstant_(maxOrder_+1),
00058 ccTerms_(maxOrder_+1),
00059 cvTerms_(maxOrder_+1),
00060 vcTerms_(maxOrder_+1),
00061 vvTerms_(maxOrder_+1),
00062 startingVectors_(maxOrder_+1),
00063 startingParities_(maxOrder_+1)
00064 {
00065 int verb = context.evalSetupVerbosity();
00066
00067 try
00068 {
00069 Tabs tabs(0);
00070 {
00071 Tabs tab;
00072 Tabs tabz;
00073 SUNDANCE_MSG1(verb,
00074 tabs << "initializing product evaluator for "
00075 << expr->toString());
00076
00077 SUNDANCE_MSG2(verb,
00078 tab << "return sparsity " << std::endl
00079 << tabz << *(this->sparsity)());
00080
00081 SUNDANCE_MSG2(verb,
00082 tab << "left sparsity " << std::endl
00083 << tabz << *(leftSparsity()) << std::endl
00084 << tabz << "right sparsity " << std::endl
00085 << tabz << *(rightSparsity()));
00086
00087 SUNDANCE_MSG3(verb,
00088 tab << "left vector index map "
00089 << leftEval()->vectorIndexMap() << std::endl
00090 << tabz << "right vector index map "
00091 << rightEval()->vectorIndexMap() << std::endl
00092 << tabz << "left constant index map "
00093 << leftEval()->constantIndexMap() << std::endl
00094 << tabz << "right constant index map "
00095 << rightEval()->constantIndexMap());
00096 }
00097 int vecResultIndex = 0;
00098 int constResultIndex = 0;
00099
00100 for (int i=0; i<this->sparsity()->numDerivs(); i++)
00101 {
00102 Tabs tab0;
00103 const MultipleDeriv& d = this->sparsity()->deriv(i);
00104
00105 SUNDANCE_MSG2(verb,
00106 tabs << std::endl
00107 << tabs << "finding rules for deriv " << d);
00108
00109 int order = d.order();
00110
00111
00112 bool resultIsConstant = this->sparsity()->state(i)==ConstantDeriv;
00113 resultIsConstant_[order].append(resultIsConstant);
00114 if (resultIsConstant)
00115 {
00116 SUNDANCE_MSG3(verb,
00117 tab0 << std::endl
00118 << tab0 << "result will be in constant index " << constResultIndex);
00119 resultIndex_[order].append(constResultIndex);
00120 addConstantIndex(i, constResultIndex);
00121 constResultIndex++;
00122 }
00123 else
00124 {
00125 SUNDANCE_MSG3(verb,
00126 tab0 << std::endl
00127 << tab0 << "result will be in constant index " << vecResultIndex);
00128 resultIndex_[order].append(vecResultIndex);
00129 addVectorIndex(i, vecResultIndex);
00130 vecResultIndex++;
00131 }
00132
00133
00134
00135
00136
00137
00138
00139 int dnLeftIndex = leftSparsity()->getIndex(d);
00140 int dnRightIndex = rightSparsity()->getIndex(d);
00141
00142
00143 bool hasVectorWorkspace = false;
00144 bool workspaceIsLeft = false;
00145 int workspaceIndex = -1;
00146 int workspaceCoeffIndex = -1;
00147 bool workspaceCoeffIsConstant = false;
00148
00149
00150 bool dnLeftIsConst = false;
00151 if (dnLeftIndex != -1)
00152 {
00153 dnLeftIsConst = leftSparsity()->state(dnLeftIndex)==ConstantDeriv;
00154 }
00155 bool dnRightIsConst = false;
00156 if (dnRightIndex != -1)
00157 {
00158 dnRightIsConst = rightSparsity()->state(dnRightIndex)==ConstantDeriv;
00159 }
00160
00161
00162 if (dnLeftIndex != -1 && !dnLeftIsConst
00163 && rightSparsity()->containsDeriv(MultipleDeriv()))
00164 {
00165
00166 hasVectorWorkspace = true;
00167 workspaceIndex = leftEval()->vectorIndexMap().get(dnLeftIndex);
00168 SUNDANCE_MSG3(verb,
00169 tab0 << "using left as workspace");
00170 workspaceIsLeft = true;
00171 int d0RightIndex = rightSparsity()->getIndex(MultipleDeriv());
00172 bool d0RightIsConst = rightSparsity()->state(d0RightIndex)==ConstantDeriv;
00173 workspaceCoeffIsConstant = d0RightIsConst;
00174 if (d0RightIsConst)
00175 {
00176 workspaceCoeffIndex
00177 = rightEval()->constantIndexMap().get(d0RightIndex);
00178 }
00179 else
00180 {
00181 workspaceCoeffIndex
00182 = rightEval()->vectorIndexMap().get(d0RightIndex);
00183 }
00184 }
00185
00186
00187 if (!hasVectorWorkspace && dnRightIndex != -1 && !dnRightIsConst
00188 && leftSparsity()->containsDeriv(MultipleDeriv()))
00189 {
00190
00191 hasVectorWorkspace = true;
00192 workspaceIndex = rightEval()->vectorIndexMap().get(dnRightIndex);
00193 workspaceIsLeft = false;
00194 SUNDANCE_MSG3(verb,
00195 tab0 << "using right as workspace");
00196 int d0LeftIndex = leftSparsity()->getIndex(MultipleDeriv());
00197 bool d0LeftIsConst = leftSparsity()->state(d0LeftIndex)==ConstantDeriv;
00198 workspaceCoeffIsConstant = d0LeftIsConst;
00199 if (d0LeftIsConst)
00200 {
00201 workspaceCoeffIndex
00202 = leftEval()->constantIndexMap().get(d0LeftIndex);
00203 }
00204 else
00205 {
00206 workspaceCoeffIndex
00207 = leftEval()->vectorIndexMap().get(d0LeftIndex);
00208 }
00209 }
00210
00211 if (hasVectorWorkspace && verb > 2)
00212 {
00213 std::string wSide = "right";
00214 MultipleDeriv wsDeriv;
00215 if (workspaceIsLeft)
00216 {
00217 wSide = "left";
00218 wsDeriv = leftSparsity()->deriv(dnLeftIndex);
00219 }
00220 else
00221 {
00222 wsDeriv = rightSparsity()->deriv(dnRightIndex);
00223 }
00224 SUNDANCE_MSG2(verb, tab0 << "has workspace vector: "
00225 << wSide << " deriv= "
00226 << wsDeriv
00227 << ", coeff index= " << workspaceCoeffIndex);
00228 }
00229 else
00230 {
00231 SUNDANCE_MSG2(verb, tab0 << "has no workspace vector");
00232 }
00233
00234 hasWorkspace_[order].append(hasVectorWorkspace);
00235 workspaceIsLeft_[order].append(workspaceIsLeft);
00236 workspaceIndex_[order].append(workspaceIndex);
00237 workspaceCoeffIndex_[order].append(workspaceCoeffIndex);
00238 workspaceCoeffIsConstant_[order].append(workspaceCoeffIsConstant);
00239
00240 ProductRulePerms perms;
00241 d.productRulePermutations(perms);
00242 SUNDANCE_MSG4(verb,
00243 tabs << "product rule permutations " << perms);
00244
00245 Array<Array<int> > ccTerms;
00246 Array<Array<int> > cvTerms;
00247 Array<Array<int> > vcTerms;
00248 Array<Array<int> > vvTerms;
00249 Array<int> startingVector;
00250 ProductParity startingParity;
00251
00252 bool hasStartingVector = false;
00253
00254 for (ProductRulePerms::const_iterator
00255 iter = perms.begin(); iter != perms.end(); iter++)
00256 {
00257 Tabs tab1;
00258 MultipleDeriv dLeft = iter->first.first();
00259 MultipleDeriv dRight = iter->first.second();
00260 int multiplicity = iter->second;
00261
00262
00263
00264 if (hasVectorWorkspace && workspaceIsLeft
00265 && dLeft.order()==order) continue;
00266 if (hasVectorWorkspace && !workspaceIsLeft
00267 && dRight.order()==order) continue;
00268
00269 if (!leftSparsity()->containsDeriv(dLeft)
00270 || !rightSparsity()->containsDeriv(dRight)) continue;
00271 int iLeft = leftSparsity()->getIndex(dLeft);
00272 int iRight = rightSparsity()->getIndex(dRight);
00273
00274 if (iLeft==-1 || iRight==-1) continue;
00275 SUNDANCE_MSG4(verb,
00276 tab1 << "left deriv=" << dLeft);
00277 SUNDANCE_MSG4(verb,
00278 tab1 << "right deriv=" << dRight);
00279
00280 bool leftIsConst = leftSparsity()->state(iLeft)==ConstantDeriv;
00281 bool rightIsConst = rightSparsity()->state(iRight)==ConstantDeriv;
00282
00283 if (leftIsConst)
00284 {
00285 Tabs tab2;
00286 SUNDANCE_MSG4(verb,
00287 tab2 << "left is const");
00288 int leftIndex = leftEval()->constantIndexMap().get(iLeft);
00289 if (rightIsConst)
00290 {
00291 SUNDANCE_MSG4(verb,
00292 tab2 << "right is const");
00293 int rightIndex = rightEval()->constantIndexMap().get(iRight);
00294 ccTerms.append(tuple(leftIndex, rightIndex, multiplicity));
00295 }
00296 else
00297 {
00298 SUNDANCE_MSG4(verb,
00299 tab2 << "right is vec");
00300 int rightIndex = rightEval()->vectorIndexMap().get(iRight);
00301 if (!hasVectorWorkspace && !hasStartingVector)
00302 {
00303 SUNDANCE_MSG4(verb,
00304 tab1 << "found c-v starting vec");
00305 startingVector = tuple(leftIndex, rightIndex,
00306 multiplicity);
00307 hasStartingVector = true;
00308 startingParity = ConstVec;
00309 }
00310 else
00311 {
00312 SUNDANCE_MSG4(verb,
00313 tab1 << "found c-v term");
00314 cvTerms.append(tuple(leftIndex, rightIndex,
00315 multiplicity));
00316 }
00317 }
00318 }
00319 else
00320 {
00321 Tabs tab2;
00322 SUNDANCE_MSG4(verb,
00323 tab2 << "left is vec");
00324 int leftIndex = leftEval()->vectorIndexMap().get(iLeft);
00325 if (rightIsConst)
00326 {
00327 SUNDANCE_MSG4(verb,
00328 tab2 << "right is const");
00329 int rightIndex = rightEval()->constantIndexMap().get(iRight);
00330 if (!hasVectorWorkspace && !hasStartingVector)
00331 {
00332 SUNDANCE_MSG4(verb,
00333 tab1 << "found v-c starting vec");
00334 startingVector = tuple(leftIndex, rightIndex,
00335 multiplicity);
00336 hasStartingVector = true;
00337 startingParity = VecConst;
00338 }
00339 else
00340 {
00341 SUNDANCE_MSG4(verb,
00342 tab1 << "found v-c term");
00343 vcTerms.append(tuple(leftIndex, rightIndex,
00344 multiplicity));
00345 }
00346 }
00347 else
00348 {
00349 SUNDANCE_MSG4(verb,
00350 tab2 << "right is vec");
00351 int rightIndex = rightEval()->vectorIndexMap().get(iRight);
00352 if (!hasVectorWorkspace && !hasStartingVector)
00353 {
00354 SUNDANCE_MSG4(verb,
00355 tab1 << "found v-v starting vec");
00356 startingVector = tuple(leftIndex, rightIndex,
00357 multiplicity);
00358 hasStartingVector = true;
00359 startingParity = VecVec;
00360 }
00361 else
00362 {
00363 SUNDANCE_MSG4(verb,
00364 tab1 << "found v-v term");
00365 vvTerms.append(tuple(leftIndex, rightIndex,
00366 multiplicity));
00367 }
00368 }
00369 }
00370 }
00371 ccTerms_[order].append(ccTerms);
00372 cvTerms_[order].append(cvTerms);
00373 vcTerms_[order].append(vcTerms);
00374 vvTerms_[order].append(vvTerms);
00375 startingVectors_[order].append(startingVector);
00376 startingParities_[order].append(startingParity);
00377
00378 if (verb > 2)
00379 {
00380 Tabs tab0;
00381 Out::os() << tab0 << "deriv " << i << " order=" << order ;
00382 if (resultIsConstant)
00383 {
00384 Out::os() << " constant result, index= ";
00385 }
00386 else
00387 {
00388 Out::os() << " vector result, index= ";
00389 }
00390 Out::os() << resultIndex_[order] << std::endl;
00391 {
00392 Tabs tab1;
00393 if (hasStartingVector)
00394 {
00395 Out::os() << tab1 << "starting vector " << startingVector << std::endl;
00396 }
00397 Out::os() << tab1 << "c-c terms " << ccTerms << std::endl;
00398 Out::os() << tab1 << "c-v terms " << cvTerms << std::endl;
00399 Out::os() << tab1 << "v-c terms " << vcTerms << std::endl;
00400 Out::os() << tab1 << "v-v terms " << vvTerms << std::endl;
00401 }
00402 }
00403 }
00404
00405 if (verb > 2)
00406 {
00407 Out::os() << tabs << "maps: " << std::endl;
00408 Out::os() << tabs << "vector index map " << vectorIndexMap() << std::endl;
00409 Out::os() << tabs << "constant index map " << constantIndexMap() << std::endl;
00410 }
00411 }
00412 catch(std::exception& e)
00413 {
00414 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
00415 "exception detected in ProductEvaluator: expr="
00416 << expr->toString() << std::endl
00417 << "exception=" << e.what());
00418 }
00419 }
00420
00421 void ProductEvaluator
00422 ::internalEval(const EvalManager& mgr,
00423 Array<double>& constantResults,
00424 Array<RCP<EvalVector> >& vectorResults) const
00425 {
00426 Tabs tabs(0);
00427
00428 SUNDANCE_MSG1(mgr.verb(),
00429 tabs << "ProductEvaluator::eval() expr="
00430 << expr()->toString());
00431
00432
00433 Array<RCP<EvalVector> > leftVectorResults;
00434 Array<RCP<EvalVector> > rightVectorResults;
00435 Array<double> leftConstantResults;
00436 Array<double> rightConstantResults;
00437 evalChildren(mgr, leftConstantResults, leftVectorResults,
00438 rightConstantResults, rightVectorResults);
00439
00440 if (mgr.verb() > 2)
00441 {
00442 Out::os() << tabs << "left operand results" << std::endl;
00443 mgr.showResults(Out::os(), leftSparsity(), leftVectorResults,
00444 leftConstantResults);
00445 Out::os() << tabs << "right operand results" << std::endl;
00446 mgr.showResults(Out::os(), rightSparsity(), rightVectorResults,
00447 rightConstantResults);
00448 }
00449
00450 constantResults.resize(this->sparsity()->numConstantDerivs());
00451 vectorResults.resize(this->sparsity()->numVectorDerivs());
00452
00453
00454
00455
00456
00457 for (int order=maxOrder_; order>=0; order--)
00458 {
00459 for (int i=0; i<resultIndex_[order].size(); i++)
00460 {
00461 double constantVal = 0.0;
00462 const Array<Array<int> >& ccTerms = ccTerms_[order][i];
00463 for (int j=0; j<ccTerms.size(); j++)
00464 {
00465
00466 constantVal += leftConstantResults[ccTerms[j][0]]
00467 * rightConstantResults[ccTerms[j][1]] * ccTerms[j][2];
00468 }
00469
00470 if (resultIsConstant_[order][i])
00471 {
00472 constantResults[resultIndex_[order][i]] = constantVal;
00473 continue;
00474 }
00475
00476
00477
00478
00479
00480
00481 RCP<EvalVector> result;
00482 if (hasWorkspace_[order][i])
00483 {
00484 int index = workspaceIndex_[order][i];
00485 int coeffIndex = workspaceCoeffIndex_[order][i];
00486 bool coeffIsConstant = workspaceCoeffIsConstant_[order][i];
00487 if (workspaceIsLeft_[order][i])
00488 {
00489 result = leftVectorResults[index];
00490 if (coeffIsConstant)
00491 {
00492 const double& coeff = rightConstantResults[coeffIndex];
00493 if (!isOne(coeff)) result->multiply_S(coeff);
00494 }
00495 else
00496 {
00497 const RCP<EvalVector>& coeff
00498 = rightVectorResults[coeffIndex];
00499 result->multiply_V(coeff.get());
00500 }
00501 }
00502 else
00503 {
00504 result = rightVectorResults[index];
00505 if (coeffIsConstant)
00506 {
00507 const double& coeff = leftConstantResults[coeffIndex];
00508 if (!isOne(coeff)) result->multiply_S(coeff);
00509 }
00510 else
00511 {
00512 const RCP<EvalVector>& coeff
00513 = leftVectorResults[coeffIndex];
00514 result->multiply_V(coeff.get());
00515 }
00516 }
00517 SUNDANCE_MSG5(mgr.verb(), tabs << "workspace = " << result->str());
00518 }
00519 else
00520 {
00521 result = mgr.popVector();
00522 const Array<int>& sv = startingVectors_[order][i];
00523 int multiplicity = sv[2];
00524 switch(startingParities_[order][i])
00525 {
00526 case VecVec:
00527 if (!isZero(constantVal))
00528 {
00529 if (!isOne(multiplicity))
00530 {
00531 result->setTo_S_add_SVV(constantVal,
00532 multiplicity,
00533 leftVectorResults[sv[0]].get(),
00534 rightVectorResults[sv[1]].get());
00535 }
00536 else
00537 {
00538 result->setTo_S_add_VV(constantVal,
00539 leftVectorResults[sv[0]].get(),
00540 rightVectorResults[sv[1]].get());
00541 }
00542 }
00543 else
00544 {
00545 if (!isOne(multiplicity))
00546 {
00547 result->setTo_SVV(multiplicity,
00548 leftVectorResults[sv[0]].get(),
00549 rightVectorResults[sv[1]].get());
00550 }
00551 else
00552 {
00553 result->setTo_VV(leftVectorResults[sv[0]].get(),
00554 rightVectorResults[sv[1]].get());
00555 }
00556 }
00557 SUNDANCE_MSG5(mgr.verb(), tabs << "init to v-v prod, m="
00558 << multiplicity << ", left="
00559 << leftVectorResults[sv[0]]->str()
00560 << ", right="
00561 << rightVectorResults[sv[1]]->str()
00562 << ", result=" << result->str());
00563 break;
00564 case VecConst:
00565 if (!isZero(constantVal))
00566 {
00567 if (!isOne(multiplicity*rightConstantResults[sv[1]]))
00568 {
00569 result->setTo_S_add_SV(constantVal,
00570 multiplicity
00571 *rightConstantResults[sv[1]],
00572 leftVectorResults[sv[0]].get());
00573 }
00574 else
00575 {
00576 result->setTo_S_add_V(constantVal,
00577 leftVectorResults[sv[0]].get());
00578 }
00579 }
00580 else
00581 {
00582 if (!isOne(multiplicity*rightConstantResults[sv[1]]))
00583 {
00584 result->setTo_SV(multiplicity
00585 *rightConstantResults[sv[1]],
00586 leftVectorResults[sv[0]].get());
00587 }
00588 else
00589 {
00590 result->setTo_V(leftVectorResults[sv[0]].get());
00591 }
00592 }
00593 SUNDANCE_MSG5(mgr.verb(), tabs << "init to v-c prod, m="
00594 << multiplicity << ", left="
00595 << leftVectorResults[sv[0]]->str()
00596 << ", right="
00597 << rightConstantResults[sv[1]]
00598 << ", result=" << result->str());
00599 break;
00600 case ConstVec:
00601 if (!isZero(constantVal))
00602 {
00603 if (!isOne(multiplicity*leftConstantResults[sv[0]]))
00604 {
00605 result->setTo_S_add_SV(constantVal,
00606 multiplicity
00607 *leftConstantResults[sv[0]],
00608 rightVectorResults[sv[1]].get());
00609 }
00610 else
00611 {
00612 result->setTo_S_add_V(constantVal,
00613 rightVectorResults[sv[1]].get());
00614 }
00615 }
00616 else
00617 {
00618 if (!isOne(multiplicity*leftConstantResults[sv[0]]))
00619 {
00620 result->setTo_SV(multiplicity
00621 *leftConstantResults[sv[0]],
00622 rightVectorResults[sv[1]].get());
00623 }
00624 else
00625 {
00626 result->setTo_V(rightVectorResults[sv[1]].get());
00627 }
00628 }
00629 SUNDANCE_MSG5(mgr.verb(), tabs << "init to c-v prod, m="
00630 << multiplicity << ", left="
00631 << leftConstantResults[sv[0]]
00632 << ", right="
00633 << rightVectorResults[sv[1]]->str()
00634 << ", result=" << result->str());
00635
00636 }
00637 SUNDANCE_MSG4(mgr.verb(), tabs << "starting value = " << result->str());
00638 }
00639 vectorResults[resultIndex_[order][i]] = result;
00640
00641 const Array<Array<int> >& cvTerms = cvTerms_[order][i];
00642 const Array<Array<int> >& vcTerms = vcTerms_[order][i];
00643 const Array<Array<int> >& vvTerms = vvTerms_[order][i];
00644
00645 for (int j=0; j<cvTerms.size(); j++)
00646 {
00647 SUNDANCE_MSG4(mgr.verb(), tabs << "adding c-v term " << cvTerms[j]);
00648
00649 double multiplicity = cvTerms[j][2];
00650 double scalar = multiplicity*leftConstantResults[cvTerms[j][0]];
00651 const EvalVector* vec = rightVectorResults[cvTerms[j][1]].get();
00652 if (!isOne(scalar))
00653 {
00654 result->add_SV(scalar, vec);
00655 }
00656 else
00657 {
00658 result->add_V(vec);
00659 }
00660 }
00661
00662 for (int j=0; j<vcTerms.size(); j++)
00663 {
00664 SUNDANCE_MSG4(mgr.verb(), tabs << "adding v-c term " << vcTerms[j]);
00665
00666 double multiplicity = vcTerms[j][2];
00667 double scalar = multiplicity*rightConstantResults[vcTerms[j][1]];
00668 const EvalVector* vec = leftVectorResults[vcTerms[j][0]].get();
00669 if (!isOne(scalar))
00670 {
00671 result->add_SV(scalar, vec);
00672 }
00673 else
00674 {
00675 result->add_V(vec);
00676 }
00677 }
00678
00679 for (int j=0; j<vvTerms.size(); j++)
00680 {
00681 SUNDANCE_MSG4(mgr.verb(), tabs << "adding v-v term " << vvTerms[j]);
00682
00683 double multiplicity = vvTerms[j][2];
00684 const EvalVector* vec1 = leftVectorResults[vvTerms[j][0]].get();
00685 const EvalVector* vec2 = rightVectorResults[vvTerms[j][1]].get();
00686 if (!isOne(multiplicity))
00687 {
00688 result->add_SVV(multiplicity, vec1, vec2);
00689 }
00690 else
00691 {
00692 result->add_VV(vec1, vec2);
00693 }
00694 }
00695
00696
00697 }
00698 }
00699
00700 if (mgr.verb() > 1)
00701 {
00702 Out::os() << tabs << "product result " << std::endl;
00703 mgr.showResults(Out::os(), this->sparsity(), vectorResults,
00704 constantResults);
00705 }
00706 }
