问题描述
S。BorağanAruoba和JesúsFernández-Villaverde编写了一套相对著名的经济学编码语言基准，here可以找到这些基准。对于C ++，有两种实现可以解决随机新古典增长模型，一种使用C样式数组，另一种使用C ++ STL数组。最近，我开始练习C ++，因为我开始在一些大型模型上工作，这种语言可能会有用（尽管我希望主要在Julia中工作），并决定使用Eigen创建此代码的变体图书馆，因为我相对不熟悉它的使用。有趣的是，我使用Eigen的代码比上面链接中的两个C ++程序都慢得多。我的假设是这是由于我的实现不佳而导致的，并且由于增加了开销，这个问题太简单了，以至于Eigen不值得使用，但是我想问一下这个假设是否正确？
我使用特征值的代码是：
#include <chrono>
#include <iostream>
#include <Eigen/Dense>
#include <cmath>
int main()
{
//---------------------------------------------------------------------------------------------------
//Timer Start
const auto time_0 = std::chrono::steady_clock::Now();
//---------------------------------------------------------------------------------------------------
//Parameter Declaration
const double aalpha = 0.33333333333; // Elasticity of output w.r.t. capital
const double bbeta = 0.95; // discount factor;
const size_t nGridProductivity = 5;
Eigen::Matrix<double,nGridProductivity,1> vProductivity;
vProductivity << 0.9792,0.9896,1.0000,1.0106,1.0212;
Eigen::Matrix<double,nGridProductivity> mTransition;
mTransition << 0.9727,0.0273,0.0000,0.0041,0.9806,0.0153,0.0082,0.9837,0.9727;
//---------------------------------------------------------------------------------------------------
// Steady State
const double capitalSteadyState = std::pow(aalpha * bbeta,1. / (1. - aalpha));
const double outputSteadyState = std::pow(capitalSteadyState,aalpha);
const double consumptionSteadyState = outputSteadyState - capitalSteadyState;
std::cout << "Output = " << outputSteadyState << ",Capital = " << capitalSteadyState << ",Consumption = " << consumptionSteadyState << "\n";
std::cout << " ";
// Capital Grid
const int nGridCapital = 17820;
Eigen::VectorXd vGridCapital(nGridCapital);
vGridCapital.setLinSpaced(nGridCapital,0.5 * capitalSteadyState,1.5 * capitalSteadyState);
//---------------------------------------------------------------------------------------------------
//required Matrices and Vectors
Eigen::MatrixXd mOutput(nGridCapital,nGridProductivity);
Eigen::MatrixXd mValueFunction (nGridCapital,nGridProductivity);
mValueFunction.setZero();//Need to initialise to zero for first iteration
Eigen::MatrixXd mValueFunctionNew(nGridCapital,nGridProductivity);
mValueFunctionNew.setZero();
Eigen::MatrixXd mPolicyFunction(nGridCapital,nGridProductivity);
mPolicyFunction.setZero();
Eigen::MatrixXd expectedValueFunction(nGridCapital,nGridProductivity);
//---------------------------------------------------------------------------------------------------
//Pre-Build Output
mOutput = (vGridCapital.array().pow(aalpha).matrix()) * vProductivity.transpose();
std::cout << "Test mOutput" << mOutput(2,2) << std::endl;
//---------------------------------------------------------------------------------------------------
//Main Iteration
const double tolerance = 0.0000001;
double maxDifference = 10.0;
std::size_t iteration = 0;
while (maxDifference > tolerance)
{
expectedValueFunction = mValueFunction * mTransition.transpose();
for (std::size_t nProductivity = 0; nProductivity < nGridProductivity; ++nProductivity)
{
// We start from prevIoUs choice (monotonicity of policy function)
std::size_t gridCapitalNextPeriod = 0;
for (std::size_t nCapital = 0; nCapital < nGridCapital; ++nCapital)
{
double valueHighSoFar = -100000.0;
double capitalChoice = vGridCapital(0);
for (std::size_t nCapitalNextPeriod = gridCapitalNextPeriod; nCapitalNextPeriod < nGridCapital; ++nCapitalNextPeriod)
{
const double consumption = mOutput(nCapital,nProductivity) - vGridCapital(nCapitalNextPeriod);
const double valueProvisional = (1. - bbeta) * std::log(consumption) + bbeta * expectedValueFunction(nCapitalNextPeriod,nProductivity);
if (valueProvisional > valueHighSoFar)
{
valueHighSoFar = valueProvisional;
capitalChoice = vGridCapital(nCapitalNextPeriod);
gridCapitalNextPeriod = nCapitalNextPeriod;
}
else
{
mValueFunctionNew(nCapital,nProductivity) = valueHighSoFar;
mPolicyFunction(nCapital,nProductivity) = capitalChoice;
// We break when we have achieved the max (note: of a monotonic function)
break;
}
mValueFunctionNew(nCapital,nProductivity) = valueHighSoFar;
mPolicyFunction(nCapital,nProductivity) = capitalChoice;
}
}
}
maxDifference = (mValueFunctionNew - mValueFunction).cwiseAbs().maxCoeff();
mValueFunction = mValueFunctionNew;
++iteration;
if ((iteration % 10 == 0) || (iteration == 1))
std::cout << "Iteration = " << iteration << ",Sup Diff = " << maxDifference << "\n";
}
std::cout << "Iteration = " << iteration << ",Sup Diff = " << maxDifference << "\n";
endl(std::cout);
std::cout << "My check = " << mPolicyFunction(999,2) << "\n";
endl(std::cout);
//---------------------------------------------------------------------------------------------------
//Timer End
const auto time_1 = std::chrono::steady_clock::Now();
const auto elapsed_seconds = std::chrono::duration_cast<std::chrono::duration<double>>(time_1 - time_0).count();
std::cout << "Elapsed time is = " << elapsed_seconds << " seconds." << std::endl;
endl(std::cout);
return 0;
}
从本质上讲，这是上面链接的代码的直接副本，其中的C数组/ STL数组被本征矩阵所取代。它需要一些整理和更多的评论，但希望它仍然可读。
在我的系统上，检查值为：
C数组版本：0.146549
STL阵列版本：0.146549
本征版本：0.146552-这很令人关注，但足够接近，以至于我可以告诉实现工作
测量的时间是：
C数组版本：0.609375
STL阵列版本：0.635014
本征版本：0.938192
程序运行之间的差异很小。
使用上面链接的github（cl /F 4000000 /o testvcpp /O2 RBC_CPP.cpp）中提到的设置，使用Visual C cl编译器对这些文件进行了编译。
解决方法
暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！
如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。
小编邮箱:dio#foxmail.com (将#修改为@）