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 "SundanceEvalVector.hpp"
00032 #include "SundanceTempStack.hpp"
00033 #include "PlayaExceptions.hpp"
00034 #include "PlayaTabs.hpp"
00035 #include "SundanceOut.hpp"
00036 #include "Teuchos_TimeMonitor.hpp"
00037 
00038 using namespace Sundance;
00039 using namespace Sundance;
00040 
00041 using namespace Teuchos;
00042 
00043 
00044 
00045 
00046 EvalVector::EvalVector(TempStack* s)
00047   : s_(s),
00048     data_(s->popVectorData()),
00049     str_()
00050 {
00051   data_->resize(s->vecSize());
00052 }
00053 
00054 
00055 EvalVector::EvalVector(TempStack* s, const RCP<Array<double> >& data,
00056                        const std::string& str)
00057   : s_(s),
00058     data_(s->popVectorData()),
00059     str_(str)
00060 {
00061   //TimeMonitor t(evalVecTimer());
00062   data_->resize(data->size());
00063   int n = data_->size();
00064 
00065   if (n > 0)
00066     {
00067       double* x = &((*data_)[0]);
00068       const double* y = &((*data)[0]);
00069       for (int i=0; i<n; i++)
00070         {
00071           x[i] = y[i];
00072         }
00073     }
00074 }
00075 
00076 EvalVector::~EvalVector()
00077 {
00078   s_->pushVectorData(data_);
00079 }
00080 
00081 
00082 void EvalVector::resize(int n)
00083 {
00084   //TimeMonitor t(evalVecTimer());
00085   data_->resize(n);
00086 }
00087 
00088 RCP<EvalVector> EvalVector::clone() const
00089 {
00090   return rcp(new EvalVector(s_, data_, str_));
00091 }
00092 
00093 void EvalVector::setToConstant(const double& alpha) 
00094 {
00095   //TimeMonitor t(evalVecTimer());
00096   int n = data_->size();
00097   if (n > 0)
00098     {
00099       double* x = &((*data_)[0]);
00100       for (int i=0; i<n; i++)
00101         {
00102           x[i] = alpha;
00103         }
00104     }
00105   if (shadowOps()) str_ = Teuchos::toString(alpha);
00106 }
00107 
00108 
00109 void EvalVector::add_SV(const double& alpha,
00110                         const EvalVector* B)
00111 {
00112   //TimeMonitor t(evalVecTimer());
00113   int n = data_->size();
00114 
00115   if (n > 0)
00116     {
00117       double* const x = start();
00118       const double* const Bx = B->start();
00119       
00120       for (int i=0; i<n; i++)
00121         {
00122           x[i] += alpha*Bx[i];
00123         }
00124 
00125       addFlops(2*n);
00126     }
00127 
00128   if (shadowOps())
00129     {
00130       if (str_ != "0") str_ = "(" + str_ + "+" 
00131         + Teuchos::toString(alpha) + "*" + B->str_ + ")";
00132       else str_ = Teuchos::toString(alpha) + "*" + B->str_;
00133     }
00134 }
00135 
00136 void EvalVector::add_S(const double& alpha)
00137 {
00138   //TimeMonitor t(evalVecTimer());
00139   int n = data_->size();
00140 
00141   if (n > 0)
00142     {
00143       double* const x = start();
00144       
00145       for (int i=0; i<n; i++)
00146         {
00147           x[i] += alpha;
00148         }
00149       addFlops(n);
00150     }
00151 
00152   if (shadowOps())
00153     {
00154       if (str_ != "0") str_ = "(" + str_ + "+" 
00155         + Teuchos::toString(alpha) + ")";
00156       else str_ = Teuchos::toString(alpha);
00157     }
00158 }
00159 
00160 
00161 void EvalVector::add_V(const EvalVector* A)
00162 {
00163   //TimeMonitor t(evalVecTimer());
00164   int n = data_->size();
00165 
00166   if (n > 0)
00167     {
00168       double* const x = start();
00169       const double* const Ax = A->start();
00170       
00171       for (int i=0; i<n; i++)
00172         {
00173           x[i] += Ax[i];
00174         }
00175       addFlops(n);
00176     }
00177 
00178   if (shadowOps())
00179     {
00180       if (str_ != "0") str_ = "(" + str_ + " + " +  A->str_ + ")";
00181       else str_ = A->str_;
00182     }
00183 }
00184 
00185 void EvalVector::add_SVV(const double& alpha,
00186                          const EvalVector* B,
00187                          const EvalVector* C)
00188 {
00189   //TimeMonitor t(evalVecTimer());
00190   int n = data_->size();
00191 
00192   if (n > 0)
00193     {
00194       double* const x = start();
00195       const double* const Bx = B->start();
00196       const double* const Cx = C->start();
00197       
00198       for (int i=0; i<n; i++)
00199         {
00200           x[i] += alpha*Bx[i]*Cx[i];
00201         }
00202       addFlops(3*n);
00203     }
00204 
00205   if (shadowOps())
00206     {
00207       if (str_ != "0") str_ = "(" + str_ + " + " 
00208         + Teuchos::toString(alpha) + "*" + B->str() + "*" + C->str() + ")";
00209       else str_ = Teuchos::toString(alpha) + "*" + B->str() + "*" + C->str();
00210     }
00211 }
00212 
00213 void EvalVector::add_VV(const EvalVector* A,
00214                         const EvalVector* B)
00215 {
00216   //TimeMonitor t(evalVecTimer());
00217   int n = data_->size();
00218 
00219   if (n > 0)
00220     {
00221       double* const x = start();
00222       const double* const Ax = A->start();
00223       const double* const Bx = B->start();
00224       
00225       for (int i=0; i<n; i++)
00226         {
00227           x[i] += Ax[i]*Bx[i];
00228         }
00229       addFlops(2*n);
00230     }
00231 
00232   if (shadowOps())
00233     {
00234       if (str_ != "0") str_ = "(" + str_ + " + " + A->str() 
00235         + "*" + B->str() + ")";
00236       else str_ =A->str() + "*" + B->str();
00237     }
00238 }
00239 
00240 
00241 void EvalVector::multiply_S_add_SV(const double& alpha, 
00242                                    const double& beta,
00243                                    const EvalVector* C)
00244 {
00245   //TimeMonitor t(evalVecTimer());
00246   int n = data_->size();
00247   
00248   if (n > 0)
00249     {
00250       double* const x = start();
00251       const double* const Cx = C->start();
00252       
00253       for (int i=0; i<n; i++)
00254         {
00255           x[i] *= alpha;
00256           x[i] += beta*Cx[i];
00257         }
00258       addFlops(3*n);
00259     }
00260 
00261   if (shadowOps())
00262     {
00263       if (str_ != "0") str_ = "(" + Teuchos::toString(alpha) + "*" + str_ + "+"
00264         + Teuchos::toString(beta) + "*" + C->str_ + ")";
00265       else str_ = Teuchos::toString(beta) + "*" + C->str_;
00266     }
00267 }
00268 
00269 
00270 void EvalVector::multiply_S_add_S(const double& alpha, 
00271                                   const double& beta)
00272 {
00273   //TimeMonitor t(evalVecTimer());
00274   int n = data_->size();
00275   
00276   if (n > 0)
00277     {
00278       double* const x = start();
00279       
00280       for (int i=0; i<n; i++)
00281         {
00282           x[i] *= alpha;
00283           x[i] += beta;
00284         }
00285       addFlops(2*n);
00286     }
00287 
00288   if (shadowOps())
00289     {
00290       if (str_ != "0") str_ = "(" + Teuchos::toString(alpha) + "*" + str_
00291         + " + " + Teuchos::toString(beta) + ")";
00292       else str_ = Teuchos::toString(beta);
00293     }
00294 }
00295 
00296 void EvalVector::multiply_V(const EvalVector* A)
00297 {
00298   //TimeMonitor t(evalVecTimer());
00299   int n = data_->size();
00300   
00301   if (n > 0)
00302     {
00303       double* const x = start();
00304       const double* const Ax = A->start();
00305       
00306       for (int i=0; i<n; i++)
00307         {
00308           x[i] *= Ax[i];
00309         }
00310       addFlops(n);
00311     }
00312 
00313   if (shadowOps())
00314     {
00315       if (str_ != "0") str_ = str() + "*" + A->str();
00316     }
00317 }
00318 
00319 void EvalVector::multiply_V_add_VVV(const EvalVector* A,
00320                                     const EvalVector* B,
00321                                     const EvalVector* C,
00322                                     const EvalVector* D)
00323 {
00324   //TimeMonitor t(evalVecTimer());
00325   int n = data_->size();
00326 
00327   if (n > 0)
00328     {
00329       double* const x = start();
00330       const double* const Ax = A->start();
00331       const double* const Bx = B->start();
00332       const double* const Cx = C->start();
00333       const double* const Dx = D->start();
00334       
00335       for (int i=0; i<n; i++)
00336         {
00337           x[i] *= Ax[i];
00338           x[i] += Bx[i]*Cx[i]*Dx[i];
00339         }
00340       addFlops(4*n);
00341     }
00342 
00343   if (shadowOps())
00344     {
00345       if (str_ != "0") str_ = "(" + str() + "*" + A->str() + " + " 
00346         + B->str() + "*" + C->str() + "*" + D->str() + ")";
00347       else str_ = B->str() + "*" + C->str() + "*" + D->str();
00348     }
00349 }
00350 
00351 void EvalVector::multiply_V_add_SVV(const EvalVector* A,
00352                                     const double& beta,
00353                                     const EvalVector* C,
00354                                     const EvalVector* D)
00355 {
00356   //TimeMonitor t(evalVecTimer());
00357   int n = data_->size();
00358 
00359   if (n > 0)
00360     {
00361       double* const x = start();
00362       const double* const Ax = A->start();
00363       const double* const Cx = C->start();
00364       const double* const Dx = D->start();
00365       
00366       for (int i=0; i<n; i++)
00367         {
00368           x[i] *= Ax[i];
00369           x[i] += beta*Cx[i]*Dx[i];
00370         }
00371       addFlops(4*n);
00372     }
00373 
00374   if (shadowOps())
00375     {
00376       if (beta != 0.0)
00377         {
00378           str_ = "(" + str() + "*" + A->str();
00379         }
00380       else
00381         {
00382           str_ = str() + "*" + A->str();
00383         }
00384       if (beta != 0.0)
00385         {
00386           str_ += " + ";
00387           if (beta != 1.0)
00388             {
00389               str_ += Teuchos::toString(beta) + "*";
00390             }
00391           str_ += C->str() + "*" + D->str();
00392         }
00393     }
00394 }
00395 
00396 void EvalVector::multiply_V_add_SV(const EvalVector* A,
00397                                    const double& beta,
00398                                    const EvalVector* C)
00399 {
00400   //TimeMonitor t(evalVecTimer());
00401   int n = data_->size();
00402 
00403   if (n > 0)
00404     {
00405       double* const x = start();
00406       const double* const Ax = A->start();
00407       const double* const Cx = C->start();
00408       
00409       for (int i=0; i<n; i++)
00410         {
00411           x[i] *= Ax[i];
00412           x[i] += beta*Cx[i];
00413         }
00414       addFlops(3*n);
00415     }
00416 
00417   if (shadowOps())
00418     {
00419       str_ = "(" + str() + "*" + A->str() + " + " + Teuchos::toString(beta)
00420         + "*" + C->str() + ")";
00421     }
00422 }
00423 
00424 void EvalVector::multiply_VV(const EvalVector* A,
00425                              const EvalVector* B)
00426 {
00427   //TimeMonitor t(evalVecTimer());
00428   int n = data_->size();
00429 
00430   if (n > 0)
00431     {
00432       double* const x = start();
00433       const double* const Ax = A->start();
00434       const double* const Bx = B->start();
00435       
00436       for (int i=0; i<n; i++)
00437         {
00438           x[i] *= Ax[i]*Bx[i];
00439         }
00440       addFlops(2*n);
00441     }
00442 
00443   if (shadowOps())
00444     {
00445       str_ = str() + "*" + A->str() + "*" + B->str();
00446     }
00447 }
00448 
00449 
00450 void EvalVector::multiply_SV(const double& alpha,
00451                              const EvalVector* B)
00452 {
00453   //TimeMonitor t(evalVecTimer());
00454   int n = data_->size();
00455 
00456   if (n > 0)
00457     {
00458       double* const x = start();
00459       
00460       const double* const Bx = B->start();
00461       
00462       for (int i=0; i<n; i++)
00463         {
00464           x[i] *= alpha*Bx[i];
00465         }
00466       addFlops(2*n);
00467     }
00468 
00469   if (shadowOps())
00470     {
00471       str_ = Teuchos::toString(alpha) + "*" + str_ + "*" + B->str();
00472     }
00473 }
00474 
00475 void EvalVector::multiply_S(const double& alpha)
00476 {
00477   //TimeMonitor t(evalVecTimer());
00478   int n = data_->size();
00479 
00480   if (n > 0)
00481     {
00482       double* const x = start();
00483       
00484       for (int i=0; i<n; i++)
00485         {
00486           x[i] *= alpha;
00487         }
00488       addFlops(n);
00489     }
00490 
00491   if (shadowOps())
00492     {
00493       str_ = Teuchos::toString(alpha) + "*" + str_;
00494     }
00495 }
00496 
00497 
00498 void EvalVector::setTo_S_add_SVV(const double& alpha,
00499                                  const double& beta,
00500                                  const EvalVector* C,
00501                                  const EvalVector* D)
00502 {
00503   //TimeMonitor t(evalVecTimer());
00504   int n = C->data_->size();
00505   resize(n);
00506 
00507   if (n > 0)
00508     {
00509       double* const x = start();
00510       const double* const Cx = C->start();
00511       const double* const Dx = D->start();
00512       
00513       for (int i=0; i<n; i++)
00514         {
00515           x[i] = alpha + beta*Cx[i]*Dx[i];
00516         }
00517       addFlops(3*n);
00518     }
00519 
00520   if (shadowOps())
00521     {
00522       str_ = "(" + Teuchos::toString(alpha) + " + " 
00523         + Teuchos::toString(beta) + "*" 
00524         + C->str() + "*" + D->str() + ")";
00525     }
00526 }
00527 
00528 void EvalVector::setTo_S_add_VV(const double& alpha,
00529                                 const EvalVector* B,
00530                                 const EvalVector* C)
00531 {
00532   //TimeMonitor t(evalVecTimer());
00533   int n = B->data_->size();
00534   resize(n);
00535 
00536   if (n > 0)
00537     {
00538       double* const x = start();
00539       const double* const Bx = B->start();
00540       const double* const Cx = C->start();
00541       
00542       for (int i=0; i<n; i++)
00543         {
00544           x[i] = alpha + Bx[i]*Cx[i];
00545         }
00546       addFlops(2*n);
00547     }
00548 
00549   if (shadowOps())
00550     {
00551       str_ = "(" + Teuchos::toString(alpha) + " + " 
00552         + B->str() + "*" + C->str() + ")";
00553     }
00554 }
00555 
00556 void EvalVector::setTo_S_add_SV(const double& alpha,
00557                                 const double& beta,
00558                                 const EvalVector* C)
00559 {
00560   //TimeMonitor t(evalVecTimer());
00561   int n = C->data_->size();
00562   resize(n);
00563 
00564   if (n > 0)
00565     {
00566       double* const x = start();
00567       const double* const Cx = C->start();
00568       
00569       for (int i=0; i<n; i++)
00570         {
00571           x[i] = alpha + beta*Cx[i];
00572         }
00573       addFlops(2*n);
00574     }
00575 
00576   if (shadowOps())
00577     {
00578       str_ = "(" + Teuchos::toString(alpha) + " + " 
00579         + Teuchos::toString(beta) + "*" 
00580         + C->str() + ")";
00581     }
00582 }
00583 
00584 
00585 
00586 void EvalVector::setTo_S_add_V(const double& alpha,
00587                                const EvalVector* B)
00588 {
00589   //TimeMonitor t(evalVecTimer());
00590   int n = B->data_->size();
00591   resize(n);
00592 
00593   if (n > 0)
00594     {
00595       double* const x = start();
00596       const double* const Bx = B->start();
00597       
00598       for (int i=0; i<n; i++)
00599         {
00600           x[i] = alpha + Bx[i];
00601         }
00602       addFlops(n);
00603     }
00604 
00605   if (shadowOps())
00606     {
00607       str_ = "(" + Teuchos::toString(alpha) + " + " 
00608         + B->str() + ")";
00609     }
00610 }
00611 
00612 
00613 void EvalVector::setTo_V(const EvalVector* A)
00614 {
00615   //TimeMonitor t(evalVecTimer());
00616   int n = A->data_->size();
00617   resize(n);
00618 
00619   if (n > 0)
00620     {
00621       double* const x = start();
00622       const double* const Ax = A->start();
00623       
00624       for (int i=0; i<n; i++)
00625         {
00626           x[i] = Ax[i];
00627         }
00628     }
00629 
00630   if (shadowOps())
00631     {
00632       str_ = A->str();
00633     }
00634 }
00635 
00636 
00637 
00638 void EvalVector::setTo_SVV(const double& alpha,
00639                            const EvalVector* B,
00640                            const EvalVector* C)
00641 {
00642   //TimeMonitor t(evalVecTimer());
00643   int n = B->data_->size();
00644   resize(n);
00645 
00646   if (n > 0)
00647     {
00648       double* const x = start();
00649       const double* const Bx = B->start();
00650       const double* const Cx = C->start();
00651 
00652       
00653       for (int i=0; i<n; i++)
00654         {
00655           x[i] = alpha*Bx[i]*Cx[i];
00656         }
00657       addFlops(2*n);
00658     }
00659 
00660   if (shadowOps())
00661     {
00662       str_ = Teuchos::toString(alpha) + "*" 
00663         + B->str() + "*" + C->str();
00664     }
00665 }
00666 
00667 void EvalVector::setTo_VV(const EvalVector* A,
00668                           const EvalVector* B)
00669 {
00670   //TimeMonitor t(evalVecTimer());
00671   int n = A->data_->size();
00672   resize(n);
00673 
00674   if (n > 0)
00675     {
00676       double* const x = start();
00677       const double* const Ax = A->start();
00678       const double* const Bx = B->start();
00679 
00680       
00681       for (int i=0; i<n; i++)
00682         {
00683           x[i] = Ax[i]*Bx[i];
00684         }
00685 
00686       addFlops(n);
00687     }
00688 
00689   if (shadowOps())
00690     {
00691       str_ = A->str() + "*" + B->str();
00692     }
00693 }
00694 
00695 void EvalVector::setTo_SV(const double& alpha,
00696                           const EvalVector* B)
00697 {
00698   //TimeMonitor t(evalVecTimer());
00699   int n = B->data_->size();
00700   resize(n);
00701 
00702   if (n > 0)
00703     {
00704       double* const x = start();
00705       const double* const Bx = B->start();
00706 
00707       
00708       for (int i=0; i<n; i++)
00709         {
00710           x[i] = alpha*Bx[i];
00711         }
00712       addFlops(n);
00713     }
00714 
00715   if (shadowOps())
00716     {
00717       str_ = Teuchos::toString(alpha) + "*" 
00718         + B->str();
00719     }
00720 }
00721 
00722 void EvalVector::applyUnaryOperator(const UnaryFunctor* func, 
00723                                     Array<RCP<EvalVector> >& opDerivs)
00724 {
00725   //TimeMonitor t(evalVecTimer());
00726   int n = data_->size();
00727 
00728   int order = opDerivs.size()-1;
00729 
00730   TEUCHOS_TEST_FOR_EXCEPTION(order < 0 || order > 2, std::runtime_error,
00731                      "illegal order=" << order << " in "
00732                      "EvalVector::applyUnaryOperator()");
00733 
00734   opDerivs[0] = s_->popVector();
00735   opDerivs[0]->resize(n);
00736   if (order > 0)
00737     {
00738       opDerivs[1] = s_->popVector();
00739       opDerivs[1]->resize(n);
00740     }
00741   if (order > 1)
00742     {
00743       opDerivs[2] = s_->popVector();
00744       opDerivs[2]->resize(n);
00745     }
00746   
00747   if (n > 0)
00748     {
00749       double* const x = start();
00750       double* const f = opDerivs[0]->start();
00751       if (order==0)
00752         {
00753           func->eval0(x, n, f);
00754         }
00755       else if (order==1)
00756         {
00757           double* const df = opDerivs[1]->start();
00758           func->eval1(x, n, f, df);
00759         }
00760       else 
00761         {
00762           double* const df = opDerivs[1]->start();
00763           double* const df2 = opDerivs[2]->start();
00764           func->eval2(x, n, f, df, df2);
00765         }
00766 
00767     }
00768   
00769   if (shadowOps())
00770     {
00771       opDerivs[0]->setString(func->name() + "(" + str() + ")");
00772       if (order > 0)
00773         {
00774           opDerivs[1]->setString(func->name() + "'(" + str() + ")");
00775         }
00776       if (order > 1)
00777         {
00778           opDerivs[2]->setString(func->name() + "\"(" + str() + ")");
00779         }
00780     }
00781 }
00782 
00783 void EvalVector::print(std::ostream& os) const 
00784 {
00785   TEUCHOS_TEST_FOR_EXCEPTION(shadowOps() && str_.size()==0, std::runtime_error, "empty eval vector result string!");
00786   os << str_;
00787 
00788   if (data_->size() > 0)
00789     {
00790       os << ", " << *data_;
00791     }
00792 }
00793 
00794 
00795 
00796 
00797 
00798 
00799 
00800
SundanceEvalVector.cpp
