stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Corrections for min, max and dummy functions in DynareBison.yy: hand_side involves hand_side arguments 

git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1420 ac1d8469-bf42-47a9-8791-bf33cf982152
time-shift
ferhat 2007-10-05 22:11:47 +00:00
parent 61dc2c742c
commit 5e8bf1ef51
7 changed files with 1803 additions and 1593 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2566
parser.src/DynareBison.cc
View File

File diff suppressed because it is too large Load Diff

10
parser.src/DynareBison.yy
View File

@ -293,8 +293,8 @@ expression : '(' expression ')'
comma_expression : expression
{ driver.add_unknown_function_arg($1); }
/* | comma_expression COMMA expression
{ driver.add_unknown_function_arg($3); }*/
| comma_expression COMMA expression
{ driver.add_unknown_function_arg($3); }
;
initval : INITVAL ';' initval_list END
@ -401,11 +401,11 @@ hand_side : '(' hand_side ')'
{ $$ = driver.add_atan($3); }
| SQRT '(' hand_side ')'
{ $$ = driver.add_sqrt($3); }
| DUMMY '(' expression ')'
| DUMMY '(' hand_side ')'
{ $$ = driver.add_dummy($3); }
| MAX '(' expression COMMA expression ')'
| MAX '(' hand_side COMMA hand_side ')'
{ $$ = driver.add_max($3 , $5); }
| MIN '(' expression COMMA expression ')'
| MIN '(' hand_side COMMA hand_side ')'
{ $$ = driver.add_min($3 , $5); }
;

22
parser.src/ModelTree.cc
View File

@ -625,7 +625,7 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl
prev_Simulation_Type=ModelBlock->Block_List[j].Simulation_Type;
}
ModelBlock_Aggregated_Count++;
cout << "ModelBlock_Aggregated_Count=" << ModelBlock_Aggregated_Count << "\n";
//cout << "ModelBlock_Aggregated_Count=" << ModelBlock_Aggregated_Count << "\n";
//For each block
j=0;
for(k0 = 0;k0 < ModelBlock_Aggregated_Count;k0++)
@ -638,7 +638,7 @@ ModelTree::writeModelEquationsCodeOrdered(const string file_name, const Model_Bl
code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
v=ModelBlock->Block_List[j].Simulation_Type;
code_file.write(reinterpret_cast<char *>(&v),sizeof(v));
cout << "FBEGINBLOCK j=" << j << " size=" << ModelBlock_Aggregated_Number[k0] << " type=" << v << "\n";
//cout << "FBEGINBLOCK j=" << j << " size=" << ModelBlock_Aggregated_Number[k0] << " type=" << v << "\n";
for(k=0; k<ModelBlock_Aggregated_Size[k0]; k++)
{
for(i=0; i < ModelBlock->Block_List[j].Size;i++)
@ -1458,28 +1458,29 @@ ModelTree::Write_Inf_To_Bin_File(const string &dynamic_basename, const string &b
}
for(j=0;j<block_triangular.ModelBlock->Block_List[num].Size;j++)
{
//int eq=block_triangular.ModelBlock->Block_List[i].Equation[j];
int eqr1=j;
int varr=block_triangular.ModelBlock->Block_List[num].Size*(block_triangular.periods
+/*block_triangular.ModelBlock->Block_List[num].Max_Lead*/block_triangular.Model_Max_Lead);
+block_triangular.Model_Max_Lead);
int k1=0;
//mDynamicModelFile << " var_in_equ_and_lag[std::make_pair(std::make_pair(" << j << ", 0)] = -1;\n";
//mDynamicModelFile << " /*periods=" << block_triangular.periods << " Size=" << block_triangular.ModelBlock->Block_List[i].Size << "*/\n";
//mDynamicModelFile << " var_in_equ_and_lag.insert(std::make_pair(std::make_pair(" << eqr1 << ", 0), " << block_triangular.ModelBlock->Block_List[i].Size*block_triangular.periods << "));\n";
//mDynamicModelFile << " equ_in_var_and_lag[std::make_pair(std::make_pair(-1, " << k1 << ")] = " << equ << ";\n";
SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1));
SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
SaveCode.write(reinterpret_cast<char *>(&k1), sizeof(k1));
SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1));
//cout << " IM_i[std::make_pair(std::make_pair(" << eqr1 << ", " << varr << "), " << k1 << ")] = " << eqr1 << ";\n";
u_count_int++;
}
//cout << "u_count_int = " << u_count_int << endl;
for(j=0;j<block_triangular.ModelBlock->Block_List[num].Size;j++)
{
//mDynamicModelFile << " index_var[" << j << "]=" << block_triangular.ModelBlock->Block_List[i].Variable[j] << ";\n";
int varr=block_triangular.ModelBlock->Block_List[num].Variable[j];
SaveCode.write(reinterpret_cast<char *>(&varr), sizeof(varr));
}
for(j=0;j<block_triangular.ModelBlock->Block_List[num].Size;j++)
{
//mDynamicModelFile << " index_var[" << j << "]=" << block_triangular.ModelBlock->Block_List[i].Variable[j] << ";\n";
int eqr1=block_triangular.ModelBlock->Block_List[num].Equation[j];
SaveCode.write(reinterpret_cast<char *>(&eqr1), sizeof(eqr1));
}
SaveCode.close();
}
@ -2413,6 +2414,7 @@ ModelTree::evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_
int j=0;
bool *IM=NULL;
int a_variable_lag=-9999;
//block_triangular.Print_IM(2);
for(first_derivatives_type::iterator it = first_derivatives.begin();
it != first_derivatives.end(); it++)
{
@ -2435,7 +2437,7 @@ ModelTree::evaluateJacobian(const eval_context_type &eval_context, jacob_map *j_
}
if (IM[eq*symbol_table.endo_nbr+var] && (fabs(val) < cutoff))
{
//cout << "the coefficient related to variable " << var << " with lag " << k1 << " in equation " << eq << " is equal to " << interprete_.u1 << " and is set to 0 in the incidence matrix\n";
//cout << "the coefficient related to variable " << var << " with lag " << k1 << " in equation " << eq << " is equal to " << val << " and is set to 0 in the incidence matrix\n";
block_triangular.unfill_IM(eq, var, k1);
i++;
}

135
parser.src/include/linbcg.hh Normal file
View File

@ -0,0 +1,135 @@
#ifndef LINBCG_HH_INCLUDED
#define LINBCG_HH_INCLUDED
#include <iostream>
#include <iomanip>
#include <cmath>
#include <map>
#include "mex.h"
using namespace std;
const int debile=10;
typedef double DP;
template <class T>
class NRVec {
private:
int nn; // size of array. upper index is nn-1
T *v;
public:
NRVec();
explicit NRVec(int n); // Zero-based array
NRVec(const T &a, int n); //initialize to constant value
NRVec(const T *a, int n); // Initialize to array
NRVec(const NRVec &rhs); // Copy constructor
NRVec & operator=(const NRVec &rhs); //assignment
NRVec & operator=(const T &a); //assign a to every element
inline T & operator[](const int i); //i'th element
inline const T & operator[](const int i) const;
inline int size() const;
~NRVec();
};
template <class T>
NRVec<T>::NRVec() : nn(0), v(0) {}
template <class T>
NRVec<T>::NRVec(int n) : nn(n), v(new T[n]) {}
template <class T>
NRVec<T>::NRVec(const T& a, int n) : nn(n), v(new T[n])
{
for(int i=0; i<n; i++)
v[i] = a;
}
template <class T>
NRVec<T>::NRVec(const T *a, int n) : nn(n), v(new T[n])
{
for(int i=0; i<n; i++)
v[i] = *a++;
}
template <class T>
NRVec<T>::NRVec(const NRVec<T> &rhs) : nn(rhs.nn), v(new T[nn])
{
for(int i=0; i<nn; i++)
v[i] = rhs[i];
}
template <class T>
NRVec<T> & NRVec<T>::operator=(const NRVec<T> &rhs)
// postcondition: normal assignment via copying has been performed;
// if vector and rhs were different sizes, vector
// has been resized to match the size of rhs
{
if (this != &rhs)
{
if (nn != rhs.nn) {
if (v != 0) delete [] (v);
nn=rhs.nn;
v= new T[nn];
}
for (int i=0; i<nn; i++)
v[i]=rhs[i];
}
return *this;
}
template <class T>
NRVec<T> & NRVec<T>::operator=(const T &a) //assign a to every element
{
for (int i=0; i<nn; i++)
v[i]=a;
return *this;
}
template <class T>
inline T & NRVec<T>::operator[](const int i) //subscripting
{
return v[i];
}
template <class T>
inline const T & NRVec<T>::operator[](const int i) const //subscripting
{
return v[i];
}
template <class T>
inline int NRVec<T>::size() const
{
return nn;
}
extern double *u;
template <class T>
NRVec<T>::~NRVec()
{
if (v != 0)
delete[] (v);
}
typedef const NRVec<int> Vec_I_INT;
typedef NRVec<int> Vec_INT, Vec_O_INT, Vec_IO_INT;
typedef const NRVec<DP> Vec_I_DP;
typedef NRVec<DP> Vec_DP, Vec_O_DP, Vec_IO_DP;
class LinBCG
{
Vec_INT *ija_p;
Vec_DP *sa_p;
public:
LinBCG();
void Init(int periods, int y_kmin, int y_kmax, int Size, std::map<std::pair<std::pair<int, int> ,int>, int> IM, int* index_vara, int* index_equa);
void asolve(Vec_I_DP &b, Vec_O_DP &x, const int itrnsp);
void atimes(Vec_I_DP &x, Vec_O_DP &r, const int itrnsp);
DP snrm(Vec_I_DP &sx, const int itol);
void sprsax(Vec_I_DP &sa, Vec_I_INT &ija, Vec_I_DP &x, Vec_O_DP &b);
void sprstx(Vec_I_DP &sa, Vec_I_INT &ija, Vec_I_DP &x, Vec_O_DP &b);
void linbcg(Vec_I_DP &b, Vec_IO_DP &x, const int itol, const DP tol,
const int itmax, int &iter, DP &err);
};
#endif

12
parser.src/include/simulate.hh
View File

@ -1,14 +1,12 @@
#ifndef SIMULATE_HH_INCLUDED
#define SIMULATE_HH_INCLUDED
#include <math>
#include <stack>
#include <set>
#include <vector>
#include <math.h>
#include <iostream>
#include <fstream>
//#include "pctimer_h.hh"
#include <time.h>
#include <string>
#include <map>
@ -20,11 +18,6 @@
#define pow_ pow
//#define pow pow1
// typedef struct IM_compact
// {
// int size, u_init, u_finish, nb_endo;
// int *u, *Var, *Equ, *Var_Index, *Equ_Index, *Var_dyn_Index;
// };
typedef struct Variable_l
{
int* Index;
@ -98,7 +91,7 @@ typedef struct NonZeroElem
std::map<std::pair<std::pair<int, int> ,int>, int> IM_i;
std::multimap<std::pair<int,int>, int> var_in_equ_and_lag_i, equ_in_var_and_lag_i;
int *index_vara, *pivota=NULL, *save_op_all=NULL, *b, *g_save_op=NULL;
int *index_vara, *index_equa, *pivota=NULL, *save_op_all=NULL, *b, *g_save_op=NULL;
int *pivot, *pivotk;
longd *pivotv, *pivotva=NULL;
bool *line_done;
@ -189,6 +182,7 @@ class SparseMatrix
NonZeroElem* NZE_Mem;
};
#define get_code_int *((int*)Code); Code+=sizeof(int);
#define get_code_double *((double*)Code); Code+=sizeof(double);
#define get_code_pdouble (double*)Code; Code+=sizeof(double);
@ -210,7 +204,7 @@ class Interpreter
{
protected :
void compute_block_time();
void simulate_a_block(int size,int type, string file_name, string bin_basename);
void simulate_a_block(int size,int type, string file_name, string bin_basename, bool Gaussian_Elimination);
longd *T;
vector<Block_contain_type> Block_Contain;
vector<Block_type> Block;

304
parser.src/linbcg.cc Normal file
View File

@ -0,0 +1,304 @@
#include "linbcg.hh"
LinBCG::LinBCG()
{
}
void LinBCG::Init(int periods, int y_kmin, int y_kmax, int Size, std::map<std::pair<std::pair<int, int> ,int>, int> IM, int* index_vara, int* index_equa)
{
std::map<std::pair<std::pair<int, int> ,int>, int>::iterator it4;
int Non_Zero_By_Lag[y_kmax+y_kmin+1];
memset(Non_Zero_By_Lag, 0, sizeof(int)*(y_kmax+y_kmin+1));
int i, j, k, nz_diag_count, eq, eqp, eq_count, var, lag, pos_u;
it4=IM.begin();
// count the non zero diagonal elements
j=0;
nz_diag_count=0;
eqp=it4->first.first.first;
mexPrintf("IM.size()=%d\n",IM.size());
for(i=0;i<Size;i++)
{
mexPrintf("index_vara[%d]=%d index_equa[%d]=%d\n",i,index_vara[i],i,index_equa[i]);
}
while (it4!=IM.end())
{
eq= it4->first.first.first;
// Assumption no Row is completly null
if(eq!=eqp)
j++;
var=it4->first.first.second;
lag=it4->first.second;
mexPrintf("eq=%d var=%d lag=%d \n",eq, var, lag);
if(!lag)
{
if(var==index_vara[j] && eq==index_equa[j])
nz_diag_count++;
else
Non_Zero_By_Lag[y_kmin+lag]++;
}
else
Non_Zero_By_Lag[y_kmin+lag]++;
it4++;
}
int Size_SparseMatrixRow=Size*periods;
mexPrintf("En lag 0 Size*periods=%d\n",Size*periods);
for(lag=-y_kmin;lag<0;lag++)
{
Size_SparseMatrixRow+=-min(-(periods+lag),0)*Non_Zero_By_Lag[y_kmin+lag];
mexPrintf("En lag %d -min(-(periods+lag),0)*Non_Zero_By_Lag[y_kmin+lag]=%d\n",lag,-min(-(periods+lag),0)*Non_Zero_By_Lag[y_kmin+lag]);
}
Size_SparseMatrixRow+=Non_Zero_By_Lag[y_kmin];
for(lag=1;lag<=y_kmax;lag++)
{
Size_SparseMatrixRow+=-min(-(periods-lag),0)*Non_Zero_By_Lag[y_kmin+lag];
mexPrintf("En lag %d -min(-(periods-lag),0)*Non_Zero_By_Lag[y_kmin+lag]=%d\n",lag,-min(-(periods-lag),0)*Non_Zero_By_Lag[y_kmin+lag]);
}
mexPrintf("Size_SparseMatrixRow=%d periods=%d\n",Size_SparseMatrixRow,periods);
//ija_p=new Vec_INT(Size_SparseMatrixRow);
//sa_p=new Vec_DP(Size_SparseMatrixRow);
Vec_INT ija_p[Size_SparseMatrixRow];
Vec_DP sa_p[Size_SparseMatrixRow];
int *cumulate_past_block;
cumulate_past_block=(int*)mxMalloc((periods+1)*sizeof(int));
memset(cumulate_past_block, 0, sizeof(int)*(periods));
mexPrintf("OK-1\n");
cumulate_past_block[0]=0;
for(j=1;j<=periods;j++)
{
for(i=-y_kmin;i<=y_kmax;i++)
{
if(j+i>0 && j+i<periods)
{
mexPrintf("cumulate_past_block[%d]+=Non_Zero_By_Lag[%d]\n",j,i+y_kmin);
cumulate_past_block[j]+=Non_Zero_By_Lag[i+y_kmin];
}
}
}
int IM_Size=IM.size();
k=Size*periods;
mexPrintf("OK00\n");
ija_p[1]=Size*periods+2;
//ija_p[Size*periods]=Size_SparseMatrixRow+1;
eq_count=1;eqp=-1;
mexPrintf("OK0\n");
it4=IM.begin();
while (it4!=IM.end())
{
eq= it4->first.first.first;
// Assumption no Row is completly null
/*if(eq!=eqp)
j++;*/
var=it4->first.first.second;
lag=it4->first.second;
pos_u=it4->second;
//mexPrintf("eq=%d var=%d lag=%d \n",eq, var, lag);
if(!lag)
{
if(var==index_vara[eq_count-1] && eq==index_equa[eq_count-1])
//nz_diag_count++;
{
//Diagonal elements
for(j=0;j<periods;j++)
{
sa_p[var+j*Size]=u[pos_u+j*IM_Size];
mexPrintf("sa_p[%d]=%f\n",var+j*Size,sa_p[var+j*Size]);
}
}
else
{
k++;
for(j=0;j<periods;j++)
{
sa_p[cumulate_past_block[j]+k]=u[pos_u+j*IM_Size];
mexPrintf("sa_p[%d]=%f\n",cumulate_past_block[j]+k,sa_p[cumulate_past_block[j]+k]);
ija_p[cumulate_past_block[j]+k]=var+1;
mexPrintf("ija_p[%d]=%d\n",cumulate_past_block[j]+k,ija_p[cumulate_past_block[j]+k]);
}
}
}
else
{
k++;
for(j=0;j<periods;j++)
if(j+lag>0 && j+lag<periods)
{
sa_p[cumulate_past_block[j]+k]=u[pos_u+j*IM_Size];
ija_p[cumulate_past_block[j]+k]=var+1;
sa_p[cumulate_past_block[j]+k]=u[pos_u+j*IM_Size];
mexPrintf("sa_p[%d]=%f\n",cumulate_past_block[j]+k,sa_p[cumulate_past_block[j]+k]);
ija_p[cumulate_past_block[j]+k]=var+1;
mexPrintf("ija_p[%d]=%d\n",cumulate_past_block[j]+k,ija_p[cumulate_past_block[j]+k]);
}
}
it4++;
if(eq!=eqp)
{
eq_count++;
for(j=0;j<periods;j++)
{
ija_p[j*Size+eq_count]=cumulate_past_block[j]+k;
mexPrintf("ija_p[%d]=%d\n",j*Size+eq_count,ija_p[j*Size+eq_count]);
}
}
eqp=eq;
}
mexPrintf("OK1\n");
for(j=1;j<Size_SparseMatrixRow;j++)
{
mexPrintf("%d ",j);
mexPrintf("ija=%d ",ija_p[j]);
mexPrintf("sa=%f\n",sa_p[j]);
}
}
void LinBCG::asolve(Vec_I_DP &b, Vec_O_DP &x, const int itrnsp)
{
int i;
int n=b.size();
for(i=0;i<n;i++) x[i]=((*sa_p)[i] != 0.0 ? b[i]/(*sa_p)[i] : b[i]);
}
void LinBCG::atimes(Vec_I_DP &x, Vec_O_DP &r, const int itrnsp)
{
if (itrnsp) sprstx(*sa_p,*ija_p,x,r);
else sprsax(*sa_p,*ija_p,x,r);
}
DP LinBCG::snrm(Vec_I_DP &sx, const int itol)
{
int i,isamax;
DP ans;
int n=sx.size();
if (itol <= 3) {
ans = 0.0;
for (i=0;i<n;i++) ans += sx[i]*sx[i];
return sqrt(ans);
} else {
isamax=0;
for (i=0;i<n;i++) {
if (fabs(sx[i]) > fabs(sx[isamax])) isamax=i;
}
return fabs(sx[isamax]);
}
}
void LinBCG::sprsax(Vec_I_DP &sa, Vec_I_INT &ija, Vec_I_DP &x, Vec_O_DP &b)
{
int i,k;
int n=x.size();
if (ija[0] != n+1)
mexPrintf("Error: \n sprsax: mismatched vector and matrix\n");
for (i=0;i<n;i++) {
b[i]=sa[i]*x[i];
for (k=ija[i];k<ija[i+1];k++) {
b[i] += sa[k]*x[ija[k]];
}
}
}
void LinBCG::sprstx(Vec_I_DP &sa, Vec_I_INT &ija, Vec_I_DP &x, Vec_O_DP &b)
{
int i,j,k;
int n=x.size();
if (ija[0] != (n+1))
mexPrintf("Error: \n mismatched vector and matrix in sprstx\n");
for (i=0;i<n;i++) b[i]=sa[i]*x[i];
for (i=0;i<n;i++) {
for (k=ija[i];k<ija[i+1];k++) {
j=ija[k];
b[j] += sa[k]*x[i];
}
}
}
void LinBCG::linbcg(Vec_I_DP &b, Vec_IO_DP &x, const int itol, const DP tol,
const int itmax, int &iter, DP &err)
{
DP ak,akden,bk,bkden=1.0,bknum,bnrm,dxnrm,xnrm,zm1nrm,znrm;
const DP EPS=1.0e-14;
int j;
int n=b.size();
Vec_DP p(n),pp(n),r(n),rr(n),z(n),zz(n);
iter=0;
atimes(x,r,0);
for (j=0;j<n;j++) {
r[j]=b[j]-r[j];
rr[j]=r[j];
}
//atimes(r,rr,0);
if (itol == 1) {
bnrm=snrm(b,itol);
asolve(r,z,0);
}
else if (itol == 2) {
asolve(b,z,0);
bnrm=snrm(z,itol);
asolve(r,z,0);
}
else if (itol == 3 || itol == 4) {
asolve(b,z,0);
bnrm=snrm(z,itol);
asolve(r,z,0);
znrm=snrm(z,itol);
} else mexPrintf("Error: \n illegal itol in linbcg\n");
cout << fixed << setprecision(6);
while (iter < itmax) {
++iter;
asolve(rr,zz,1);
for (bknum=0.0,j=0;j<n;j++) bknum += z[j]*rr[j];
if (iter == 1) {
for (j=0;j<n;j++) {
p[j]=z[j];
pp[j]=zz[j];
}
} else {
bk=bknum/bkden;
for (j=0;j<n;j++) {
p[j]=bk*p[j]+z[j];
pp[j]=bk*pp[j]+zz[j];
}
}
bkden=bknum;
atimes(p,z,0);
for (akden=0.0,j=0;j<n;j++) akden += z[j]*pp[j];
ak=bknum/akden;
atimes(pp,zz,1);
for (j=0;j<n;j++) {
x[j] += ak*p[j];
r[j] -= ak*z[j];
rr[j] -= ak*zz[j];
}
asolve(r,z,0);
if (itol == 1)
err=snrm(r,itol)/bnrm;
else if (itol == 2)
err=snrm(z,itol)/bnrm;
else if (itol == 3 || itol == 4) {
zm1nrm=znrm;
znrm=snrm(z,itol);
if (fabs(zm1nrm-znrm) > EPS*znrm) {
dxnrm=fabs(ak)*snrm(p,itol);
err=znrm/fabs(zm1nrm-znrm)*dxnrm;
} else {
err=znrm/bnrm;
continue;
}
xnrm=snrm(x,itol);
if (err <= 0.5*xnrm) err /= xnrm;
else {
err=znrm/bnrm;
continue;
}
}
cout << "iter=" << setw(4) << iter+1 << setw(12) << err << endl;
if (err <= tol) break;
}
}

347
parser.src/simulate.cc
View File

@ -4,6 +4,7 @@
// use NO_COMPILER option in MODEL command //
////////////////////////////////////////////////////////////////////////
#include "simulate.hh"
#include "linbcg.hh"
longd
@ -29,6 +30,34 @@ max(int a, int b)
return(b);
}
int
min(int a, int b)
{
if (a<b)
return(a);
else
return(b);
}
double
max(double a, double b)
{
if (a>b)
return(a);
else
return(b);
}
double
min(double a, double b)
{
if (a<b)
return(a);
else
return(b);
}
#ifdef INDIRECT_SIMULATE
void
read_file_table_u(t_table_u **table_u, t_table_u **F_table_u, t_table_u **i_table_u, t_table_u **F_i_table_u, int *nb_table_u, bool i_to_do, bool shifting, int *nb_add_u_count)
@ -1415,6 +1444,7 @@ void Read_SparseMatrix(std::string file_name, int Size, int periods, int y_kmin,
for (j=0;j<Size;j++)
{
SaveCode.read(reinterpret_cast<char *>(&index_vara[j]), sizeof(*index_vara));
//mexPrintf("index_vara[%d]=%d\n",j,index_vara[j]);
}
for (i=1;i<periods+y_kmin+y_kmax;i++)
{
@ -1429,6 +1459,12 @@ void Read_SparseMatrix(std::string file_name, int Size, int periods, int y_kmin,
#endif
}
}
index_equa=(int*)mxMalloc(Size*sizeof(int));
for(j=0;j<Size;j++)
{
SaveCode.read(reinterpret_cast<char *>(&index_equa[j]), sizeof(*index_equa));
//mexPrintf("index_equa[%d]=%d\n",j,index_equa[j]);
}
}
@ -3510,283 +3546,17 @@ simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, in
return(0);
}
/*
void SparseMatrixRow::Init(int periods, int y_kmin, int y_kmax, int Size, std::map<std::pair<std::pair<int, int> ,int>, int> IM, double* sa, int* ija)
{
std::map<std::pair<std::pair<int, int> ,int>, int>::iterator it4;
int Non_Zero_Elem=IM.size();
it4=IM.begin();
// count the non zero diagonal elements
while (it4!=IM.end())
{
eq= it4->first.first.first;
var=it4->first.first.second;
lag=it4->first.second;
for (t=0;t<periods;t++)
{
#ifdef PRINT_OUT
mexPrintf("t=%d\n",t);
#endif
int ti_y_kmin=-min( t , y_kmin);
int ti_y_kmax= min( periods-(t+1), y_kmax);
it4=IM.begin();
eq=-1;
while (it4!=IM.end())
{
var=it4->first.first.second;
if (eq!=it4->first.first.first+Size*t)
tmp_b=0;
eq=it4->first.first.first+Size*t;
#ifdef PRINT_OUT
mexPrintf("eq=%d, var=%d",eq,var);
#endif
if (var<(periods+y_kmax)*Size)
{
lag=it4->first.second;
#ifdef PRINT_OUT
mexPrintf(", lag =%d, ti_y_kmin=%d, ti_y_kmax=%d ", lag, ti_y_kmin, ti_y_kmax);
#endif
if (lag<=ti_y_kmax && lag>=ti_y_kmin)
{
//mexPrintf("u_index=%d, eq=%d, var=%d, lag=%d ",it4->second+u_count_init*t, eq, var, lag);
var+=Size*t;
//mexPrintf("u_index=%d, eq=%d, var=%d, lag=%d ",it4->second+u_count_init*t, eq, var, lag);
NbNZRow[eq]++;
NbNZCol[var]++;
#ifdef NEW_ALLOC
first=mxMalloc_NZE();
#else
first=(NonZeroElem*)mxMalloc(sizeof(NonZeroElem));
#endif
first->NZE_C_N=NULL;
first->NZE_R_N=NULL;
first->u_index=it4->second+u_count_init*t;
first->r_index=eq;
first->c_index=var;
first->lag_index=lag;
if (FNZE_R[eq]==NULL)
{
FNZE_R[eq]=first;
}
if (FNZE_C[var]==NULL)
FNZE_C[var]=first;
if (temp_NZE_R[eq]!=NULL)
temp_NZE_R[eq]->NZE_R_N=first;
if (temp_NZE_C[var]!=NULL)
temp_NZE_C[var]->NZE_C_N=first;
temp_NZE_R[eq]=first;
temp_NZE_C[var]=first;
#ifdef PRINT_OUT
mexPrintf("=> ");
#endif
}
else
{
#ifdef PRINT_OUT
mexPrintf("nn ");
mexPrintf("tmp_b+=u[%d]*y[index_var[%d]]\n",it4->second+u_count_init*t,var+Size*(y_kmin+t));
mexPrintf("tmp_b+=u[%d](%f)*y[%d(%d)](%f)",it4->second+u_count_init*t,u[it4->second+u_count_init*t], index_vara[var+Size*(y_kmin+t)],var+Size*(y_kmin+t),y[index_vara[var+Size*(y_kmin+t)]]);
#endif
tmp_b+=u[it4->second+u_count_init*t]*y[index_vara[var+Size*(y_kmin+t)]];
}
}
else
{
#ifdef PRINT_OUT
mexPrintf("");
#endif
b[eq]=it4->second+u_count_init*t;
u[b[eq]]+=tmp_b;
#ifdef PRINT_OUT
mexPrintf("=> b");
#endif
}
#ifdef PRINT_OUT
mexPrintf(" u[%d] = %e\n",it4->second+u_count_init*t,double(u[it4->second+u_count_init*t]));
#endif
it4++;
}
}
int t,i, eq, var, lag;
longd tmp_b=0.0;
NonZeroElem* first;
#ifdef MEM_ALLOC_CHK
mexPrintf("pivot=(int*)mxMalloc(%d*sizeof(int))\n",Size*periods);
#endif
pivot=(int*)mxMalloc(Size*periods*sizeof(int));
#ifdef MEM_ALLOC_CHK
mexPrintf("pivota=(int*)mxMalloc(%d*sizeof(int))\n",Size*periods);
#endif
pivotk=(int*)mxMalloc(Size*periods*sizeof(int));
#ifdef MEM_ALLOC_CHK
mexPrintf("pivotv=(longd*)mxMalloc(%d*sizeof(longd))\n",Size*periods);
#endif
pivotv=(longd*)mxMalloc(Size*periods*sizeof(longd));
#ifdef MEM_ALLOC_CHK
mexPrintf("pivotva=(longd*)mxMalloc(%d*sizeof(longd))\n",Size*periods);
#endif
pivotva=(longd*)mxMalloc(Size*periods*sizeof(longd));
#ifdef MEM_ALLOC_CHK
mexPrintf("b=(int*)mxMalloc(%d*sizeof(int))\n",Size*periods);
#endif
b=(int*)mxMalloc(Size*periods*sizeof(int));
#ifdef MEM_ALLOC_CHK
mexPrintf("line_done=(bool*)mxMalloc(%d*sizeof(bool))\n",Size*periods);
#endif
line_done=(bool*)mxMalloc(Size*periods*sizeof(bool));
memset(line_done, 0, periods*Size*sizeof(*line_done));
CHUNK_BLCK_SIZE=u_count;
print_err=true;
swp_f=false;
g_save_op=NULL;
g_nop_all=0;
#ifdef PRINT_OUT
mexPrintf("sizeof(NonZeroElem)=%d sizeof(NonZeroElem*)=%d\n",sizeof(NonZeroElem),sizeof(NonZeroElem*));
#endif
i=(periods+y_kmax+1)*Size*sizeof(NonZeroElem*);
#ifdef MEM_ALLOC_CHK
mexPrintf("FNZE_R=(NonZeroElem**)mxMalloc(%d)\n",i);
#endif
FNZE_R=(NonZeroElem**)mxMalloc(i);
#ifdef MEM_ALLOC_CHK
mexPrintf("FNZE_C=(NonZeroElem**)mxMalloc(%d)\n",i);
#endif
FNZE_C=(NonZeroElem**)mxMalloc(i);
memset(FNZE_R, 0, i);
memset(FNZE_C, 0, i);
#ifdef MEM_ALLOC_CHK
mexPrintf("temp_NZE_R=(NonZeroElem**)(%d)\n",i);
#endif
NonZeroElem** temp_NZE_R=(NonZeroElem**)mxMalloc(i);
#ifdef MEM_ALLOC_CHK
mexPrintf("temp_NZE_R=(NonZeroElem**)(%d)\n",i);
#endif
NonZeroElem** temp_NZE_C=(NonZeroElem**)mxMalloc(i);
memset(temp_NZE_R, 0, i);
memset(temp_NZE_C, 0, i);
i=(periods+y_kmax+1)*Size*sizeof(int);
#ifdef MEM_ALLOC_CHK
mexPrintf("NbNZRow=(int*)mxMalloc(%d)\n",i);
#endif
NbNZRow=(int*)mxMalloc(i);
#ifdef MEM_ALLOC_CHK
mexPrintf("NbNZCol=(int*)mxMalloc(%d)\n",i);
#endif
NbNZCol=(int*)mxMalloc(i);
#ifdef MEM_ALLOC_CHK
mexPrintf("ok\n");
#endif
memset(NbNZRow, 0, i);
memset(NbNZCol, 0, i);
i=periods*Size*sizeof(*b);
memset(b,0,i);
#ifdef PRINT_OUT
mexPrintf("Now looping\n");
#endif
for (t=0;t<periods;t++)
{
#ifdef PRINT_OUT
mexPrintf("t=%d\n",t);
#endif
int ti_y_kmin=-min( t , y_kmin);
int ti_y_kmax= min( periods-(t+1), y_kmax);
it4=IM.begin();
eq=-1;
while (it4!=IM.end())
{
var=it4->first.first.second;
if (eq!=it4->first.first.first+Size*t)
tmp_b=0;
eq=it4->first.first.first+Size*t;
#ifdef PRINT_OUT
mexPrintf("eq=%d, var=%d",eq,var);
#endif
if (var<(periods+y_kmax)*Size)
{
lag=it4->first.second;
#ifdef PRINT_OUT
mexPrintf(", lag =%d, ti_y_kmin=%d, ti_y_kmax=%d ", lag, ti_y_kmin, ti_y_kmax);
#endif
if (lag<=ti_y_kmax && lag>=ti_y_kmin)
{
//mexPrintf("u_index=%d, eq=%d, var=%d, lag=%d ",it4->second+u_count_init*t, eq, var, lag);
var+=Size*t;
//mexPrintf("u_index=%d, eq=%d, var=%d, lag=%d ",it4->second+u_count_init*t, eq, var, lag);
NbNZRow[eq]++;
NbNZCol[var]++;
#ifdef NEW_ALLOC
first=mxMalloc_NZE();
#else
first=(NonZeroElem*)mxMalloc(sizeof(NonZeroElem));
#endif
first->NZE_C_N=NULL;
first->NZE_R_N=NULL;
first->u_index=it4->second+u_count_init*t;
first->r_index=eq;
first->c_index=var;
first->lag_index=lag;
if (FNZE_R[eq]==NULL)
{
FNZE_R[eq]=first;
}
if (FNZE_C[var]==NULL)
FNZE_C[var]=first;
if (temp_NZE_R[eq]!=NULL)
temp_NZE_R[eq]->NZE_R_N=first;
if (temp_NZE_C[var]!=NULL)
temp_NZE_C[var]->NZE_C_N=first;
temp_NZE_R[eq]=first;
temp_NZE_C[var]=first;
#ifdef PRINT_OUT
mexPrintf("=> ");
#endif
}
else
{
#ifdef PRINT_OUT
mexPrintf("nn ");
mexPrintf("tmp_b+=u[%d]*y[index_var[%d]]\n",it4->second+u_count_init*t,var+Size*(y_kmin+t));
mexPrintf("tmp_b+=u[%d](%f)*y[%d(%d)](%f)",it4->second+u_count_init*t,u[it4->second+u_count_init*t], index_vara[var+Size*(y_kmin+t)],var+Size*(y_kmin+t),y[index_vara[var+Size*(y_kmin+t)]]);
#endif
tmp_b+=u[it4->second+u_count_init*t]*y[index_vara[var+Size*(y_kmin+t)]];
}
}
else
{
#ifdef PRINT_OUT
mexPrintf("");
#endif
b[eq]=it4->second+u_count_init*t;
u[b[eq]]+=tmp_b;
#ifdef PRINT_OUT
mexPrintf("=> b");
#endif
}
#ifdef PRINT_OUT
mexPrintf(" u[%d] = %e\n",it4->second+u_count_init*t,double(u[it4->second+u_count_init*t]));
#endif
it4++;
}
}
mxFree(temp_NZE_R);
mxFree(temp_NZE_C);
}
*/
/*class EvalException
{
};*/
Interpreter::Interpreter()
{
GaussSeidel=true;
//GaussSeidel=true;
}
void
@ -4187,7 +3957,7 @@ Interpreter::compute_block_time() /*throw(EvalException)*/
}
void
Interpreter::simulate_a_block(int size,int type, string file_name, string bin_basename)
Interpreter::simulate_a_block(int size,int type, string file_name, string bin_basename, bool Gaussian_Elimination)
{
char *begining;
int i;
@ -4199,8 +3969,10 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
int giter;
int u_count_int;
double *y_save;
LinBCG linbcg;
//mexPrintf("%d\n",debile);
GaussSeidel=false;
//GaussSeidel=false;
//slowc_save=slowc/2;
switch (type)
{
@ -4546,7 +4318,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
mexPrintf("u_count=%d\n",u_count);
#endif
begining=get_code_pointer;
GaussSeidel=false;
//GaussSeidel=false;
giter=0;
//mexPrintf("GaussSeidel=%d\n",GaussSeidel);
if (!is_linear)
@ -4571,7 +4343,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
if (isnan(res1)||isinf(res1))
{
memcpy(y, y_save, y_size*sizeof(longd)*(periods+y_kmax+y_kmin));
GaussSeidel=false;
//GaussSeidel=false;
break;
}
for (i=0; i< size; i++)
@ -4585,32 +4357,27 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
max_res=fabs(rr);
res2+=rr*rr;
res1+=fabs(rr);
if (GaussSeidel && giter)
/*if (GaussSeidel && giter)
{
//mexPrintf("y[%d]-=r[%d]/u[%d]\n",Block_Contain[i].Variable,i,Block_Contain[i].Own_Derivative,);
y[Per_y_+Block_Contain[i].Variable]-=r[i]/u[Per_u_+Block_Contain[i].Own_Derivative];
//mexPrintf("y[%d]-=r[%d](%f)/u[%d](%f)=%f\n",Block_Contain[i].Variable,i,r[i],Block_Contain[i].Own_Derivative,u[Per_u_+Block_Contain[i].Own_Derivative], y[Per_y_+Block_Contain[i].Variable]);
}
}*/
}
}
cvg=(max_res<solve_tolf);
if (GaussSeidel)
if(Gaussian_Elimination)
simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg);
else
{
mexPrintf("GaussSeidel res1=%f res2=%f max_res=%f\n",res1,res2,max_res);
if (giter && res1>res1a)
{
memcpy(y, y_save, nb_endo*sizeof(longd)*(periods+y_kmax+y_kmin));
GaussSeidel=false;
}
giter++;
res1a=res1;
}
if (!GaussSeidel)
{
simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg);
//GaussSeidel=true;
iter++;
int ITMAX=75;
double TOL=1e-9;
int ITOL=1;
linbcg.Init(periods, y_kmin, y_kmax, size, IM_i, index_vara, index_equa);
/*Vec_DP b(NP),x(NP);
linbcg.linbcg(b,x,ITOL,TOL,ITMAX,iter,err);*/
}
iter++;
}
if (!cvg)
{
@ -4641,7 +4408,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
mxFree(u);
mxFree(index_vara);
memset(direction,0,size_of_direction);
GaussSeidel=false;
//GaussSeidel=false;
break;
default:
mexPrintf("Unknow type =%d\n",type);
@ -4732,7 +4499,7 @@ Interpreter::compute_blocks(string file_name, string bin_basename)
//mexPrintf("Block_Contain[%d].Own_Derivative=%d\n",i,lBlock_Contain.Own_Derivative);
Block_Contain.push_back(lBlock_Contain);
}
simulate_a_block(lBlock.size,lBlock.type, file_name, bin_basename);
simulate_a_block(lBlock.size,lBlock.type, file_name, bin_basename,/*false*/true);
break;
case FEND :
//mexPrintf("FEND\n");