Commit b7f96a47 authored by george's avatar george
Browse files

- adding overloaded multInvRight to PLUFact,

- inlining log likelihood and 
- refactoring class constructors from kalman loop


git-svn-id: https://www.dynare.org/svn/dynare/trunk@2842 ac1d8469-bf42-47a9-8791-bf33cf982152
parent a52d38c8
......@@ -688,6 +688,7 @@ BasicKalmanTask::filterNonDiffuse(const Vector&a,const GeneralMatrix&P,
GeneralMatrix PtLttrans(m,m);
GeneralMatrix Mt(m,n);
GeneralMatrix Kt(m,n);
GeneralMatrix KtTransTmp(n,m); // perm space for temp Kt trans
GeneralMatrix QtRttrans(rcols,Rt.numRows());
bool isTunit=0;// Tt->isUnit();
......@@ -697,8 +698,11 @@ BasicKalmanTask::filterNonDiffuse(const Vector&a,const GeneralMatrix&P,
const double neg_alpha=-1.0;
const double omega=0.0;
Vector vt(n);
Vector Finvv(n);
Vector atsave(m);//((const Vector&)at);
int p;
int t= 1;
double vFinvv,ll;
int nonsteady=1;
for(;t<=data.numCols()&&nonsteady;++t)
{
......@@ -738,7 +742,8 @@ BasicKalmanTask::filterNonDiffuse(const Vector&a,const GeneralMatrix&P,
// GeneralMatrix Kt(Tt,Mt);
BLAS_dgemm("N", "N", &m, &n, &m, &alpha, Tt.base(), &m,
Mt.base(), &m, &omega, Kt.base(), &m);
Ftinv.multInvRight(Kt);
// Ftinv.multInvRight(Kt);
Ftinv.multInvRight(Kt, KtTransTmp); // using perm space for temp KT trans
// Kt.multInvRight(Ft);
/*****************
......@@ -755,8 +760,14 @@ BasicKalmanTask::filterNonDiffuse(const Vector&a,const GeneralMatrix&P,
/*****************
Here we calc likelihood and store results.
*****************/
double ll= calcStepLogLik(Ftinv,vt);
// fres.set(t,Ftinv,vt,Lt,at,Pt,ll);
// double ll= calcStepLogLik(Ftinv,vt);
p= Ftinv.numRows();
Finvv=vt;
Ftinv.multInvLeft(Finvv);
vFinvv= vt.dot(Finvv);
ll=-0.5*(p*log(2*M_PI)+Ftinv.getLogDeterminant()+vFinvv);
// fres.set(t,Ftinv,vt,Lt,at,Pt,ll);
(*vll)[t-1]=ll;
if (t>start) loglik+=ll;
......
......@@ -196,13 +196,26 @@ void PLUFact::multInvRight(GeneralMatrix&a)const
{
GeneralMatrix atrans(a,"trans");
TS_RAISE_IF(rows!=atrans.numRows(),
"Wrong dimension of the matrix in PLUFact::multInvLeft");
"Wrong dimension of the matrix in PLUFact::multInvRight");
PL_dgetrs("T",atrans.getData().base(),atrans.getLD(),atrans.numCols());
for(int i= 0;i<a.numRows();i++)
for(int j= 0;j<a.numCols();j++)
a.get(i,j)= atrans.get(j,i);
}
// pass also a temporary GM space for atrans to avoid matrix construction:
void PLUFact::multInvRight(GeneralMatrix&a, GeneralMatrix&atrans)const
{
TS_RAISE_IF(rows!=atrans.numRows(),
"Wrong dimension of the matrix in PLUFact::multInvRight");
for(int i= 0;i<a.numRows();i++)
for(int j= 0;j<a.numCols();j++)
atrans.get(j,i)= a.get(i,j);
PL_dgetrs("T",atrans.getData().base(),atrans.getLD(),atrans.numCols());
for(int i= 0;i<a.numRows();i++)
for(int j= 0;j<a.numCols();j++)
a.get(i,j)= atrans.get(j,i);
}
;
void PLUFact::multInvLeft(Vector&a)const
......
......@@ -67,6 +67,8 @@ class PLUFact{
const PLUFact& operator = (const GeneralMatrix&m);
void multInvLeft(GeneralMatrix&a)const;
void multInvRight(GeneralMatrix&a)const;
// pass temporary GM space for atrans to avoid matrix construction:
void multInvRight(GeneralMatrix&a, GeneralMatrix&atrans)const;
void multInvLeft(Vector&a)const;
void multInvRight(Vector&a)const;
bool isRegular()const
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment