Sundance 2.1

00001 /* @HEADER@ */
00002 // ************************************************************************
00003 // 
00004 //                              Sundance
00005 //                 Copyright (2005) Sandia Corporation
00006 // 
00007 // Copyright (year first published) Sandia Corporation.  Under the terms 
00008 // of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government 
00009 // retains certain rights in this software.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //                                                                                 
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA                                                                                
00025 // Questions? Contact Kevin Long (krlong@sandia.gov), 
00026 // Sandia National Laboratories, Livermore, California, USA
00027 // 
00028 // ************************************************************************
00029 /* @HEADER@ */
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) /* case where the left operand is zero */
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) /* case where the right operand is zero */
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 /* both are nonzero */
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   //  TimeMonitor timer(evalTimer());
00206   Tabs tabs(0);
00207 
00208   SUNDANCE_MSG1(mgr.verb(),
00209                tabs << "SumEvaluator::eval() expr=" << expr()->toString());
00210 
00211   /* evaluate the children */
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   /* Do constant terms with left=0 */
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   /* Do constant terms with right=0 */
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   /* Do vector terms with left=0 */
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   /* Do vector terms with right=0 */
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   /** Do constant-constant terms */
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   /** Do constant-vector sums */
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   /* Do vector-constant sums */
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   /* Do vector-vector sums */
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
SundanceSumEvaluator.cpp
