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 "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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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