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 "SundanceSumEvaluator.hpp"
00032 #include "SundanceEvalManager.hpp"
00033 #include "SundanceSumExpr.hpp"
00034 #include "SundanceProductExpr.hpp"
00035
00036 #include "PlayaTabs.hpp"
00037 #include "SundanceOut.hpp"
00038
00039 using namespace Sundance;
00040 using namespace Sundance;
00041 using namespace Sundance;
00042 using namespace Teuchos;
00043
00044
00045
00046
00047 SumEvaluator::SumEvaluator(const SumExpr* se,
00048 const EvalContext& context)
00049 : BinaryEvaluator<SumExpr>(se, context),
00050 sign_(se->sign()),
00051 singleRightConstant_(),
00052 singleRightVector_(),
00053 singleLeftConstant_(),
00054 singleLeftVector_(),
00055 ccSums_(),
00056 cvSums_(),
00057 vcSums_(),
00058 vvSums_()
00059 {
00060 Tabs tabs(0);
00061 int verb = context.evalSetupVerbosity();
00062
00063 if (verb > 1)
00064 {
00065 Out::os() << tabs << "initializing sum evaluator for "
00066 << se->toString() << std::endl;
00067 }
00068
00069 int constantCounter = 0;
00070 int vectorCounter = 0;
00071
00072 if (verb > 2)
00073 {
00074 Out::os() << std::endl << tabs << "return sparsity ";
00075 this->sparsity()->print(Out::os());
00076 Out::os() << std::endl << tabs << "left sparsity ";
00077 leftSparsity()->print(Out::os());
00078 Out::os() << std::endl << tabs << "right sparsity " ;
00079 rightSparsity()->print(Out::os());
00080 Out::os() << std::endl;
00081
00082 Out::os() << tabs << "left vector index map " << leftEval()->vectorIndexMap() << std::endl;
00083 Out::os() << tabs << "right vector index map " << rightEval()->vectorIndexMap() << std::endl;
00084
00085 Out::os() << tabs << "left constant index map " << leftEval()->constantIndexMap() << std::endl;
00086 Out::os() << tabs << "right constant index map " << rightEval()->constantIndexMap() << std::endl;
00087 }
00088
00089 for (int i=0; i<this->sparsity()->numDerivs(); i++)
00090 {
00091 const MultipleDeriv& d = this->sparsity()->deriv(i);
00092 TEUCHOS_TEST_FOR_EXCEPTION(!leftSparsity()->containsDeriv(d)
00093 && !rightSparsity()->containsDeriv(d),
00094 std::logic_error,
00095 "deriv " << d.toString()
00096 << " was not found in either left or right operand "
00097 "of expr " << expr()->toString());
00098
00099 int iLeft = -1;
00100 int iRight = -1;
00101
00102 if (leftSparsity()->containsDeriv(d))
00103 {
00104 iLeft = leftSparsity()->getIndex(d);
00105 }
00106
00107 if (rightSparsity()->containsDeriv(d))
00108 {
00109 iRight = rightSparsity()->getIndex(d);
00110 }
00111
00112 SUNDANCE_MSG2(verb, tabs << "deriv " << d);
00113
00114 if (iLeft == -1)
00115 {
00116 SUNDANCE_MSG3(verb, tabs << "left operand is zero ");
00117 if (rightSparsity()->state(iRight)==ConstantDeriv)
00118 {
00119 int rc = rightEval()->constantIndexMap().get(iRight);
00120 singleRightConstant_.append(tuple(constantCounter, rc));
00121 addConstantIndex(i, constantCounter++);
00122 SUNDANCE_MSG4(verb, tabs << "single right constant "
00123 << singleRightConstant_[singleRightConstant_.size()-1]);
00124 }
00125 else
00126 {
00127 int rv = rightEval()->vectorIndexMap().get(iRight);
00128 singleRightVector_.append(tuple(vectorCounter, rv));
00129 addVectorIndex(i, vectorCounter++);
00130 SUNDANCE_MSG4(verb, tabs << "single right vector "
00131 << singleRightVector_[singleRightVector_.size()-1]);
00132 }
00133 }
00134 else if (iRight == -1)
00135 {
00136 SUNDANCE_MSG3(verb, tabs << "right operand is zero ");
00137 if (leftSparsity()->state(iLeft)==ConstantDeriv)
00138 {
00139 int lc = leftEval()->constantIndexMap().get(iLeft);
00140 singleLeftConstant_.append(tuple(constantCounter, lc));
00141 addConstantIndex(i, constantCounter++);
00142 SUNDANCE_MSG4(verb, tabs << "single left constant "
00143 << singleLeftConstant_[singleLeftConstant_.size()-1]);
00144 }
00145 else
00146 {
00147 int lv = leftEval()->vectorIndexMap().get(iLeft);
00148 singleLeftVector_.append(tuple(vectorCounter, lv));
00149 addVectorIndex(i, vectorCounter++);
00150 SUNDANCE_MSG4(verb, tabs << "single left vector "
00151 << singleLeftVector_[singleLeftVector_.size()-1]);
00152 }
00153 }
00154 else
00155 {
00156 SUNDANCE_MSG4(verb, tabs << "both operands are nonzero ");
00157 bool leftIsConstant = leftSparsity()->state(iLeft)==ConstantDeriv;
00158 bool rightIsConstant = rightSparsity()->state(iRight)==ConstantDeriv;
00159
00160 if (leftIsConstant && rightIsConstant)
00161 {
00162 SUNDANCE_MSG4(verb, tabs << "both operands are constant");
00163 int lc = leftEval()->constantIndexMap().get(iLeft);
00164 int rc = rightEval()->constantIndexMap().get(iRight);
00165 ccSums_.append(tuple(constantCounter, lc, rc));
00166 addConstantIndex(i, constantCounter++);
00167 SUNDANCE_MSG4(verb, tabs << "c-c sum " << ccSums_[ccSums_.size()-1]);
00168 }
00169 else if (leftIsConstant)
00170 {
00171 SUNDANCE_MSG4(verb, tabs << "left operand is constant");
00172 int lc = leftEval()->constantIndexMap().get(iLeft);
00173 int rv = rightEval()->vectorIndexMap().get(iRight);
00174 cvSums_.append(tuple(vectorCounter, lc, rv));
00175 addVectorIndex(i, vectorCounter++);
00176 SUNDANCE_MSG4(verb, tabs << "c-v sum " << cvSums_[cvSums_.size()-1]);
00177 }
00178 else if (rightIsConstant)
00179 {
00180 SUNDANCE_MSG4(verb, tabs << "right operand is constant");
00181 int lv = leftEval()->vectorIndexMap().get(iLeft);
00182 int rc = rightEval()->constantIndexMap().get(iRight);
00183 vcSums_.append(tuple(vectorCounter, lv, rc));
00184 addVectorIndex(i, vectorCounter++);
00185 SUNDANCE_MSG4(verb, tabs << "v-c sum " << vcSums_[vcSums_.size()-1]);
00186 }
00187 else
00188 {
00189 SUNDANCE_MSG4(verb, tabs << "both operands are vectors");
00190 int lv = leftEval()->vectorIndexMap().get(iLeft);
00191 int rv = rightEval()->vectorIndexMap().get(iRight);
00192 vvSums_.append(tuple(vectorCounter, lv, rv));
00193 addVectorIndex(i, vectorCounter++);
00194 SUNDANCE_MSG4(verb, tabs << "v-v sum " << vvSums_[vvSums_.size()-1]);
00195 }
00196 }
00197 }
00198 }
00199
00200 void SumEvaluator
00201 ::internalEval(const EvalManager& mgr,
00202 Array<double>& constantResults,
00203 Array<RCP<EvalVector> >& vectorResults) const
00204 {
00205
00206 Tabs tabs(0);
00207
00208 SUNDANCE_MSG1(mgr.verb(),
00209 tabs << "SumEvaluator::eval() expr=" << expr()->toString());
00210
00211
00212 Array<RCP<EvalVector> > leftVectorResults;
00213 Array<RCP<EvalVector> > rightVectorResults;
00214 Array<double> leftConstResults;
00215 Array<double> rightConstResults;
00216 evalChildren(mgr, leftConstResults, leftVectorResults,
00217 rightConstResults, rightVectorResults);
00218
00219 if (verb() > 2)
00220 {
00221 Out::os() << tabs << "left operand " << std::endl;
00222 mgr.showResults(Out::os(), leftSparsity(), leftVectorResults,
00223 leftConstResults);
00224 Out::os() << tabs << "right operand " << std::endl;
00225 mgr.showResults(Out::os(), rightSparsity(), rightVectorResults,
00226 rightConstResults);
00227 }
00228
00229 constantResults.resize(this->sparsity()->numConstantDerivs());
00230 vectorResults.resize(this->sparsity()->numVectorDerivs());
00231
00232
00233 for (int i=0; i<singleRightConstant_.size(); i++)
00234 {
00235 Tabs tab1;
00236 constantResults[singleRightConstant_[i][0]]
00237 = sign_*rightConstResults[singleRightConstant_[i][1]];
00238 SUNDANCE_MSG2(mgr.verb(), tab1 << "sum for "
00239 << constantResultDeriv(singleRightConstant_[i][0])
00240 << ": L=0, R="
00241 << sign_*rightConstResults[singleRightConstant_[i][1]]);
00242 }
00243
00244
00245 for (int i=0; i<singleLeftConstant_.size(); i++)
00246 {
00247 Tabs tab1;
00248 constantResults[singleLeftConstant_[i][0]]
00249 = leftConstResults[singleLeftConstant_[i][1]];
00250 SUNDANCE_MSG2(mgr.verb(), tab1 << "sum for "
00251 << constantResultDeriv(singleLeftConstant_[i][0])
00252 << ": L="
00253 << leftConstResults[singleLeftConstant_[i][1]]
00254 << " R=0");
00255 }
00256
00257
00258 for (int i=0; i<singleRightVector_.size(); i++)
00259 {
00260 Tabs tab1;
00261 if (sign_ < 0.0) rightVectorResults[singleRightVector_[i][1]]->multiply_S(sign_);
00262 vectorResults[singleRightVector_[i][0]]
00263 = rightVectorResults[singleRightVector_[i][1]];
00264 SUNDANCE_MSG2(mgr.verb(), tab1 << "sum for "
00265 << vectorResultDeriv(singleRightVector_[i][0])
00266 << ": L=0, R="
00267 << sign_ << "*" <<
00268 rightVectorResults[singleRightVector_[i][1]]->str());
00269 }
00270
00271
00272 for (int i=0; i<singleLeftVector_.size(); i++)
00273 {
00274 Tabs tab1;
00275 vectorResults[singleLeftVector_[i][0]]
00276 = leftVectorResults[singleLeftVector_[i][1]];
00277 SUNDANCE_MSG2(mgr.verb(), tab1 << "sum for "
00278 << vectorResultDeriv(singleLeftVector_[i][0])
00279 << ": L="
00280 << leftVectorResults[singleLeftVector_[i][1]]->str()
00281 << " R=0");
00282 }
00283
00284
00285 for (int i=0; i<ccSums_.size(); i++)
00286 {
00287 Tabs tab1;
00288 constantResults[ccSums_[i][0]]
00289 = leftConstResults[ccSums_[i][1]]
00290 + sign_*rightConstResults[ccSums_[i][2]];
00291 SUNDANCE_MSG2(mgr.verb(), tab1 << "c-c sum for "
00292 << constantResultDeriv(ccSums_[i][0])
00293 << ": L=" << leftConstResults[ccSums_[i][1]]
00294 << " R=" << sign_*rightConstResults[ccSums_[i][2]]);
00295 }
00296
00297
00298 for (int i=0; i<cvSums_.size(); i++)
00299 {
00300 Tabs tab1;
00301 RCP<EvalVector>& v = rightVectorResults[cvSums_[i][2]];
00302 SUNDANCE_MSG2(mgr.verb(), tab1 << "doing c-v sum for "
00303 << vectorResultDeriv(cvSums_[i][0])
00304 << ": L=" << leftConstResults[cvSums_[i][1]]
00305 << " R=" << sign_ << "*"
00306 << rightVectorResults[cvSums_[i][2]]->str());
00307 if (isOne(sign_))
00308 {
00309 v->add_S(leftConstResults[cvSums_[i][1]]);
00310 }
00311 else
00312 {
00313 v->multiply_S_add_S(sign_, leftConstResults[cvSums_[i][1]]);
00314 }
00315 vectorResults[cvSums_[i][0]] = v;
00316 }
00317
00318
00319 for (int i=0; i<vcSums_.size(); i++)
00320 {
00321 Tabs tab1;
00322 RCP<EvalVector>& v = leftVectorResults[vcSums_[i][1]] ;
00323 SUNDANCE_MSG2(mgr.verb(), tab1 << "doing v-c sum for "
00324 << vectorResultDeriv(vcSums_[i][0])
00325 << ": L=" << leftVectorResults[vcSums_[i][1]]->str()
00326 << " R="
00327 << sign_*rightConstResults[vcSums_[i][2]]);
00328 v->add_S(sign_*rightConstResults[vcSums_[i][2]]);
00329 vectorResults[vcSums_[i][0]] = v;
00330 }
00331
00332
00333 for (int i=0; i<vvSums_.size(); i++)
00334 {
00335 Tabs tab1;
00336 RCP<EvalVector>& v = leftVectorResults[vvSums_[i][1]];
00337 SUNDANCE_MSG2(mgr.verb(), tab1 << "doing v-v sum for "
00338 << vectorResultDeriv(vvSums_[i][0])
00339 << ": L="
00340 << leftVectorResults[vvSums_[i][1]]->str()
00341 << " R=" << sign_ << "*"
00342 << rightVectorResults[vvSums_[i][2]]->str());
00343 if (isOne(sign_))
00344 {
00345 v->add_V(rightVectorResults[vvSums_[i][2]].get());
00346 }
00347 else
00348 {
00349 v->add_SV(sign_, rightVectorResults[vvSums_[i][2]].get());
00350 }
00351 vectorResults[vvSums_[i][0]] = v;
00352 }
00353
00354 if (mgr.verb() > 1)
00355 {
00356 Out::os() << tabs << "sum result " << std::endl;
00357 mgr.showResults(Out::os(), this->sparsity(), vectorResults,
00358 constantResults);
00359 }
00360 }
00361
00362