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 "SundanceUserDefOpCommonEvaluator.hpp"
00032 #include "SundanceEvalManager.hpp"
00033 #include "SundanceEvalVector.hpp"
00034 #include "SundanceUserDefOp.hpp"
00035 #include "SundanceUserDefOpElement.hpp"
00036 #include "SundanceSparsitySuperset.hpp"
00037 #include "PlayaTabs.hpp"
00038 #include "SundanceOut.hpp"
00039
00040 using namespace Sundance;
00041 using namespace Sundance;
00042
00043 using namespace Sundance;
00044 using namespace Teuchos;
00045
00046
00047 UserDefOpCommonEvaluator
00048 ::UserDefOpCommonEvaluator(const UserDefFunctor* functor,
00049 const UserDefOpElement* expr,
00050 const EvalContext& context)
00051 : maxOrder_(0),
00052 argValueIndex_(functor->domainDim(), -1),
00053 argValueIsConstant_(functor->domainDim()),
00054 constArgDerivCache_(functor->rangeDim()),
00055 varArgDerivCache_(functor->rangeDim()),
00056 cacheIsValid_(false),
00057 functor_(functor)
00058 {
00059
00060 for (int i=0; i<functor->domainDim(); i++)
00061 {
00062 const SparsitySuperset* sArg = expr->evaluatableChild(i)->sparsitySuperset(context).get();
00063 int numConst=0;
00064 int numVec=0;
00065 for (int j=0; j<sArg->numDerivs(); j++)
00066 {
00067 if (sArg->deriv(j).order() == 0)
00068 {
00069 if (sArg->state(j)==VectorDeriv)
00070 {
00071 argValueIndex_[i] = numVec;
00072 }
00073 else
00074 {
00075 argValueIndex_[i] = numConst;
00076 }
00077 break;
00078 }
00079 if (sArg->state(j) == VectorDeriv)
00080 {
00081 numVec++;
00082 }
00083 else
00084 {
00085 numConst++;
00086 }
00087 }
00088
00089 TEUCHOS_TEST_FOR_EXCEPTION(argValueIndex_[i]==-1, std::runtime_error,
00090 "no zeroth derivative found for argument #" << i
00091 << " of " << expr->toString());
00092 }
00093 }
00094
00095
00096
00097
00098 void UserDefOpCommonEvaluator
00099 ::evalAllComponents(const EvalManager& mgr,
00100 const Array<RCP<Array<double> > >& constArgVals,
00101 const Array<RCP<Array<RCP<EvalVector> > > >& vArgVals) const
00102 {
00103 Tabs tab0;
00104 int numPoints = EvalManager::stack().vecSize();
00105 SUNDANCE_MSG3(mgr.verb(), tab0 << "UDOpCommonEval::evalAllComponents()");
00106 SUNDANCE_MSG3(mgr.verb(), tab0 << "num points = " << numPoints);
00107 SUNDANCE_MSG2(mgr.verb(), tab0 << "max diff order = " << maxOrder_);
00108
00109 TEUCHOS_TEST_FOR_EXCEPTION(numPoints==0, std::logic_error,
00110 "Empty vector detected in evalArgDerivs()");
00111
00112
00113
00114 Array<RCP<EvalVector> > argVals(argValueIndex_.size());
00115 Array<double*> argPtrs(argValueIndex_.size());
00116 for (int q=0; q<argValueIndex_.size(); q++)
00117 {
00118 Tabs tab1;
00119 if (argValueIsConstant_[q])
00120 {
00121 argVals[q] = mgr.popVector();
00122 double* ptr = argVals[q]->start();
00123 double c = (*(constArgVals[q]))[argValueIndex_[q]];
00124 for (int p=0; p<numPoints; p++)
00125 {
00126 ptr[p] = c;
00127 }
00128 argPtrs[q] = ptr;
00129 }
00130 else
00131 {
00132 argVals[q] = (*(vArgVals[q]))[argValueIndex_[q]];
00133 argPtrs[q] = argVals[q]->start();
00134 }
00135 SUNDANCE_MSG3(mgr.verb(), tab1 << "argument #" << q << " is:");
00136 Tabs tab2;
00137 SUNDANCE_MSG3(mgr.verb(), tab2 << argVals[q]->str());
00138 }
00139
00140
00141 TEUCHOS_TEST_FOR_EXCEPTION(maxOrder_ > 2, std::runtime_error,
00142 "Differentiation order " << maxOrder_ << ">2 not supported "
00143 "for user-defined operators");
00144 int rangeDim = functor_->rangeDim();
00145 int domainDim = functor_->domainDim();
00146 int nTotal = 1;
00147 int numFirst = domainDim;
00148 int numSecond = domainDim*(domainDim+1)/2;
00149 if (maxOrder_ > 0) nTotal += numFirst;
00150 if (maxOrder_ > 1) nTotal += numSecond;
00151 int numResultVecs = nTotal * rangeDim;
00152
00153
00154
00155
00156 Array<double*> resultVecs(numResultVecs);
00157
00158
00159 for (int i=0; i<rangeDim; i++)
00160 {
00161 varArgDerivCache_[i].resize(nTotal);
00162 varArgDerivCache_[i][0] = mgr.popVector();
00163 varArgDerivCache_[i][0]->resize(numPoints);
00164 varArgDerivCache_[i][0]->setString(functor_->name(i));
00165 int d0Pos = i;
00166 SUNDANCE_MSG3(mgr.verb(), "zeroth deriv of elem #" << i << " is at " << d0Pos);
00167 resultVecs[d0Pos] = varArgDerivCache_[i][0]->start();
00168 if (maxOrder_ > 0)
00169 {
00170 int ptr = 0;
00171 for (int j=0; j<domainDim; j++)
00172 {
00173 varArgDerivCache_[i][j+1] = mgr.popVector();
00174 varArgDerivCache_[i][j+1]->resize(numPoints);
00175 int d1Pos = rangeDim + domainDim*i + j;
00176 SUNDANCE_MSG3(mgr.verb(), "first deriv (" << j << ") of elem #" << i
00177 << " is at " << d1Pos);
00178 resultVecs[d1Pos]
00179 = varArgDerivCache_[i][j+1]->start();
00180 varArgDerivCache_[i][j+1]->setString("D[" + functor_->name(i)
00181 + ", "
00182 + argVals[j]->str() + "]");
00183 if (maxOrder_ > 1)
00184 {
00185 for (int k=0; k<=j; k++, ptr++)
00186 {
00187 int m = (1 + numFirst);
00188 varArgDerivCache_[i][m+ptr] = mgr.popVector();
00189 varArgDerivCache_[i][m+ptr]->resize(numPoints);
00190 varArgDerivCache_[i][m+ptr]->setString("D["
00191 + functor_->name(i)
00192 + ", {"
00193 + argVals[j]->str()
00194 + ", "
00195 + argVals[k]->str()
00196 + "}]");
00197 int d2Pos = rangeDim + domainDim*rangeDim
00198 + i*numSecond + ptr;
00199 SUNDANCE_MSG3(mgr.verb(), "second deriv (" << j << ", " << k <<
00200 ") of elem #" << i << " is at " << d2Pos);
00201 resultVecs[d2Pos]
00202 = varArgDerivCache_[i][m+ptr]->start();
00203 }
00204 }
00205 }
00206 }
00207 }
00208
00209
00210
00211
00212 const double** in = const_cast<const double**>(&(argPtrs[0]));
00213 double** out = &(resultVecs[0]);
00214
00215 functor_->evaluationCallback(numPoints, maxOrder_, in, out);
00216
00217
00218
00219 markCacheAsValid();
00220 }
00221
00222