diff --git a/mex/sources/build_matlab.m b/mex/sources/build_matlab.m
index d1f063c19355c5179c39458d3dcb5d65bfdaa104..1b02660b0f910a324c0dcceee1170f03d07b209a 100644
--- a/mex/sources/build_matlab.m
+++ b/mex/sources/build_matlab.m
@@ -99,4 +99,4 @@ eval([ COMPILE_COMMAND ' -DMATLAB -Igensylv/cc ' ...
'gensylv/cc/Vector.cpp ' ...
BLAS_PATH ' ' LAPACK_PATH ]);
disp('Compiling simulate...')
-eval([ COMPILE_COMMAND ' -Isimulate -I../../preprocessor/include simulate/simulate.cc simulate/Interpreter.cc simulate/Mem_Mngr.cc simulate/SparseMatrix.cc simulate/linbcg.cc' ]);
+eval([ COMPILE_COMMAND ' -Isimulate -I../../preprocessor/include simulate/simulate.cc simulate/Interpreter.cc simulate/Mem_Mngr.cc simulate/SparseMatrix.cc' ]);
diff --git a/mex/sources/build_octave.m b/mex/sources/build_octave.m
index 673edb2f73e834542b6f96e7715e0c9d7bd51d74..8ce220ad4020801b32c6400a6360909c61825285 100644
--- a/mex/sources/build_octave.m
+++ b/mex/sources/build_octave.m
@@ -52,4 +52,4 @@ eval([ COMPILE_COMMAND ' -DMATLAB -Igensylv/cc ' ...
'gensylv/cc/Vector.cpp ' ...
'-o ../octave/gensylv.mex']);
disp('Compiling simulate...')
-eval([ COMPILE_COMMAND ' -Isimulate -I../../preprocessor/include simulate/simulate.cc simulate/Interpreter.cc simulate/Mem_Mngr.cc simulate/SparseMatrix.cc simulate/linbcg.cc -o ../octave/simulate.mex']);
+eval([ COMPILE_COMMAND ' -Isimulate -I../../preprocessor/include simulate/simulate.cc simulate/Interpreter.cc simulate/Mem_Mngr.cc simulate/SparseMatrix.cc -o ../octave/simulate.mex']);
diff --git a/mex/sources/simulate/Interpreter.cc b/mex/sources/simulate/Interpreter.cc
index 53ab86daa47bdd9a636af469787f05b1e4c5c336..739b41c2e02787f2eb1004ce50c5c7f46253716f 100644
--- a/mex/sources/simulate/Interpreter.cc
+++ b/mex/sources/simulate/Interpreter.cc
@@ -579,9 +579,11 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
int giter;
int u_count_int;
double *y_save;
+#ifdef LINBCG
LinBCG linbcg;
Mat_DP a;
Vec_INT indx;
+#endif
//SparseMatrix sparse_matrix;
int nb_endo, u_count_init;
@@ -960,6 +962,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
begining=get_code_pointer;
if(!Gaussian_Elimination)
{
+#ifdef LINBCG
it_=y_kmin;
Per_u_=0;
Per_y_=it_*y_size;
@@ -967,6 +970,7 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
compute_block_time();
linbcg.Initialize(filename, res1, res2, max_res, slowc, ya, direction, iter);
linbcg.Preconditioner(periods, y_kmin, y_kmax, size, IM_i, index_vara, index_equa, y_size, y, true, 0, a, indx);
+#endif
}
//GaussSeidel=false;
giter=0;
@@ -1044,8 +1048,10 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
}
else
{
+#ifdef LINBCG
linbcg.Initialize(filename, res1, res2, max_res, slowc, ya, direction, iter);
linbcg.SolveLinear(periods, y_kmin, y_kmax, size, IM_i, index_vara, index_equa,y_size,y, true, cvg, a, indx);
+#endif
}
iter++;
}
@@ -1077,8 +1083,10 @@ Interpreter::simulate_a_block(int size,int type, string file_name, string bin_ba
simulate_NG1(Block_Count, symbol_table_endo_nbr, it_, y_kmin, y_kmax, size, periods, true, cvg, iter);
else
{
+#ifdef LINBCG
linbcg.Initialize(filename, res1, res2, max_res, slowc, ya, direction, iter);
linbcg.SolveLinear(periods, y_kmin, y_kmax, size, IM_i, index_vara, index_equa, y_size, y, true, cvg, a, indx);
+#endif
}
}
#ifdef DEBUGC
diff --git a/mex/sources/simulate/Interpreter.hh b/mex/sources/simulate/Interpreter.hh
index a03cb41ed99a7fc62c60ea03e0a5365f253b387a..25cd73cb13b79c8e8b7f45a21a683c7a9a6764c3 100644
--- a/mex/sources/simulate/Interpreter.hh
+++ b/mex/sources/simulate/Interpreter.hh
@@ -29,7 +29,9 @@
//#include "ExprNode.hh"
//#include "Mem_Mngr.hh"
#include "SparseMatrix.hh"
-#include "linbcg.hh"
+#ifdef LINBCG
+# include "linbcg.hh"
+#endif
#include "mex.h"
//#define DEBUGC
diff --git a/mex/sources/simulate/linbcg.cc b/mex/sources/simulate/linbcg.cc
deleted file mode 100644
index dba8e786b3e959669774706e7c67d86ccf884a48..0000000000000000000000000000000000000000
--- a/mex/sources/simulate/linbcg.cc
+++ /dev/null
@@ -1,1498 +0,0 @@
-/*
- * Copyright (C) 2007-2008 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see <http://www.gnu.org/licenses/>.
- */
-
-#include <cstring>
-
-#include "linbcg.hh"
-//#define DEBUG
-
-LinBCG::LinBCG()
-{
-}
-
-
-void LinBCG::Preconditioner(int periods, int y_kmin, int y_kmax, int Size, std::map<std::pair<std::pair<int, int> ,int>, int> IM_i, int* index_vara, int* index_equa, int y_size, double *y, bool print_it, int type, Mat_O_DP &a, Vec_O_INT &indx)
-{
- std::map<std::pair<std::pair<int, int>, int> , int>::iterator it4_i;
- std::map<std::pair< int, int >, double>::iterator it4;
- int i, j, k, eq=0, eqp, var, lag, pos_u;
- std::map<std::pair< int, int >, double> IM;
- LinBCG U_linbcg, L_linbcg;
- Vec_INT ijL, ijU;
- Vec_DP sL, sU;
- double d;
-
- IM.clear();
- it4_i=IM_i.begin();
- int u_size=IM_i.size();
- //mexPrintf("Preconditioner\n");
- double u_pos_u;
- while(it4_i!=IM_i.end())
- {
- eq= it4_i->first.first.first;
- var=it4_i->first.first.second;
- lag=it4_i->first.second;
- pos_u=it4_i->second;
- if(pos_u>=Size) // to eliminate the residual contain in the first (Size) elements of IM
- {
- //mexPrintf("IM_i[%d, (%d, %d)]=%d val=%f ",eq, var, lag, pos_u, u[pos_u]);
- IM[make_pair(eq,var-lag*Size)]+=u[pos_u/*+(periods-1)*u_size*/];
- //mexPrintf(" IM[%d, %d]=%f\n",eq,var-lag*Size, IM[make_pair(eq,var-lag*Size)]);
- }
- it4_i++;
- }
- int Size_SparseMatrixRow=IM.size();
- //mexPrintf("Size_SparseMatrixRow+1=%d\n",Size_SparseMatrixRow+1);
- int ija_d[Size_SparseMatrixRow+1];
- DP sa_d[Size_SparseMatrixRow+1];
- k=Size;
- ija_d[0]=k+1;
- j=0;
- eqp=-99;
- it4=IM.begin();
- while(it4!=IM.end())
- {
- eq=it4->first.first;
- var=it4->first.second;
- u_pos_u=it4->second;
- //mexPrintf("IM[eq=%d, var=%d]=%f \n",eq, var, u_pos_u);
- if(eq!=eqp)
- ija_d[eq+j*Size]=k+1;
- if(eq==var)
- sa_d[var+j*Size]+=u_pos_u;//u[pos_u+j*u_size];
- else
- {
- k++;
- sa_d[k]+=u_pos_u;//u[pos_u+j*u_size];
- ija_d[k]=var+j*Size;
- }
- eqp=eq;
- it4++;
- }
- eq=Size;//Size*periods;
- ija_d[eq]=k+1;
- sa_d[eq]=0;
- ija_p.Format(ija_d,Size_SparseMatrixRow+1);
- sa_p.Format(sa_d,Size_SparseMatrixRow+1);
- /*sprsprt();
- sprs_sprt();*/
- //mexPrintf("ok0\n");
- a.Format(Size,Size);
- sprsout(a);
- indx.Format(Size);
- //mexPrintf("ok1\n");
- /*Mat_DP aa(0.0,Size*periods,Size*periods);
- for(k=0;k<periods;k++)
- for(i=0;i<Size;i++)
- for(j=0;j<Size;j++)
- aa[i+k*Size][j+k*Size]=a[i][j];
- a=aa;*/
- /*mexPrintf("a before=\n");
- a.Print();*/
- sprsludcmp(a,indx, d);
- /*mexPrintf("a after=\n");
- a.Print();*/
- //mexPrintf("end of preconditioner\n");
- /*sprs_LU(sL, ijL, sU, ijU);
-
- U_linbcg.sprsin(sU, ijU);
-
- mexPrintf("U\n");
- U_linbcg.sprs_sprt();
- U_linbcg.sprsprt();
-
- mexPrintf("ijL[0]=%d\n",ijL[0]);
- L_linbcg.sprsin(sL, ijL);
- mexPrintf("L\n");
- L_linbcg.sprs_sprt();
- L_linbcg.sprsprt();*/
-}
-
-
-void LinBCG::SolveLinear(int periods, int y_kmin, int y_kmax, int Size, std::map<std::pair<std::pair<int, int> ,int>, int> IM_i, int* index_vara, int* index_equa, int y_size, double *y, bool print_it, bool cvg, Mat_DP &a, Vec_IO_INT &indx)
-{
- // Solve iteratively the systme A*x=b
- // A is sparse and rewritten using sa_p and ija_p
- mexPrintf("SolveLinear\n");
- mexEvalString("drawnow;");
- std::map<std::pair<std::pair<int, int>, int> , int>::iterator it4_i;
- std::map<std::pair<int, std::pair<int,int> >, int>::iterator it4;
- int i, j, k, eq=0, eqp, var, lag, pos_u;
- std::map<std::pair< int, std::pair<int,int> >, int> IM;
- clock_t t1 = clock();
- IM.clear();
- it4_i=IM_i.begin();
- int u_size=IM_i.size();
- int Size_SparseMatrixRow=0;
- Vec_DP b(Size*periods), x(Size*periods);
- int ITMAX=10*Size*periods;
- double TOL=1e-9;
- int ITOL=1;
- int sub_iter;
- DP err;
- //Vec_DP Err(Size*periods);
-
- if (iter>0)
- mexPrintf("Sim : %f ms\n",(1000.0*(double(clock())-double(time00)))/double(CLOCKS_PER_SEC));
- if (isnan(res1) || isinf(res1))
- {
- if (slowc_save<1e-8)
- {
- mexPrintf("Dynare cannot improve the simulation\n");
- mexEvalString("drawnow;");
- filename+=" stopped";
- mexEvalString("st=fclose('all');clear all;");
- mexErrMsgTxt(filename.c_str());
- }
- slowc_save/=2;
- mexPrintf("Error: Simulation diverging, trying to correct it using slowc=%f\n",slowc_save);
- for (i=0;i<y_size*(periods+y_kmin);i++)
- y[i]=ya[i]+slowc_save*direction[i];
- iter--;
- return;
- }
- mexPrintf("-----------------------------------\n");
- mexPrintf(" Simulate iteration� %d \n",iter+1);
- mexPrintf(" max. error=%.10e \n",double(max_res));
- mexPrintf(" sqr. error=%.10e \n",double(res2));
- mexPrintf(" abs. error=%.10e \n",double(res1));
- mexPrintf("-----------------------------------\n");
- mexPrintf("cvg=%d\n\n",cvg);
- mexEvalString("drawnow;");
- if (cvg)
- return;
- for(j=0;j<periods;j++)
- for(i=0;i<Size;i++)
- {
- //mexPrintf("b[%d]=-u[%d]=%f\n",i+j*Size,i+j*u_size,-u[i+j*u_size]);
- b[i+j*Size]=-u[i+j*u_size];
- }
- for(i=0;i<Size;i++)
- for(j=y_kmin;j<periods+y_kmin;j++)
- {
- //mexPrintf("x[%d]=y[%d]=%f\n",i+(j-y_kmin)*Size, index_vara[i]+j*y_size, y[index_vara[i]+j*y_size]);
- x[i+(j-y_kmin)*Size]=y[index_vara[i]+j*y_size];
- }
- Size_SparseMatrixRow=periods*(u_size-Size);
- while(it4_i!=IM_i.end())
- {
- eq= it4_i->first.first.first;
- var=it4_i->first.first.second;
- lag=it4_i->first.second;
- pos_u=it4_i->second;
- if(pos_u>=Size) // to eliminate the residual contain in the first (Size) elements of IM
- {
- //mexPrintf("IM[%d, (%d, %d)]=%d\n",eq, var, lag, pos_u);
- IM[make_pair(eq,make_pair(var,lag))]=pos_u;
- for(j=0;j<=min(y_kmax+y_kmin,periods-1);j++)
- {
- //mexPrintf("j=%d\n",j);
- if(j-y_kmin<=0)
- k=j;
- else
- k=periods+j-y_kmin-y_kmax-1;
- if(lag+j<0 || lag+j>min(periods-1,y_kmax+y_kmin))
- {
- //mexPrintf("eliminate b[eq+(k)*Size=%d]=u[%d]*y[%d]\n",eq+(k)*Size,pos_u+(k)*u_size,var+(k+1)*y_size);
- b[eq+(k)*Size]-=u[pos_u+k*u_size]*y[index_vara[var-lag*Size]+(k+lag+1)*y_size];
- Size_SparseMatrixRow--;
- }
-
- }
- }
- it4_i++;
- }
- //mexPrintf("Size_SparseMatrixRow+1=%d\n",Size_SparseMatrixRow+1);
- int ija_d[Size_SparseMatrixRow+1];
- DP sa_d[Size_SparseMatrixRow+1];
- k=Size*periods;
- ija_d[0]=k+1;
- for(j=0;j<periods;j++)
- {
- eqp=-99;
- it4=IM.begin();
- while(it4!=IM.end())
- {
- eq=it4->first.first;
- var=it4->first.second.first;
- lag=it4->first.second.second;
- pos_u=it4->second;
- /*if(j==1)
- mexPrintf("IM[eq=%d, var=%d, lag=%d]=%d val=%f \n",eq, var, lag, pos_u,u[pos_u+j*u_size]);*/
- if(lag+j>=0 && lag+j<periods)
- {
- if(eq!=eqp)
- {
- ija_d[eq+j*Size]=k+1;
- //mexPrintf("ija_d[%d](1)=%d\n",eq+j*Size,k+1);
- }
- if(eq==var)
- {
-#ifdef DEBUG
- sa_d[var+j*Size]=double(pos_u+j*u_size);
-#else
- sa_d[var+j*Size]=u[pos_u+j*u_size];
-#endif
- //mexPrintf("sa_d[%d]=u[%d]=%f\n",var+j*Size,pos_u+j*u_size,u[pos_u+j*u_size]);
- }
- else
- {
- if(j+lag>=0 && j+lag<periods)
- {
- k++;
-#ifdef DEBUG
- sa_d[k]=double(pos_u+j*u_size);
-#else
- sa_d[k]=u[pos_u+j*u_size];
-#endif
- //mexPrintf("sa_d[%d]=u[%d]=%f\n",k,pos_u+j*u_size,u[pos_u+j*u_size]);
- ija_d[k]=var+j*Size;
- //mexPrintf("ija_d[%d](3)=%d\n",k,var+j*Size);
- }
- }
- eqp=eq;
- }
- it4++;
- }
- }
- eq=Size*periods;
- //mexPrintf("eq=%d\n",eq);
- ija_d[eq]=k+1;
- sa_d[eq]=0;
- /*for(j=0;j<Size_SparseMatrixRow+1;j++)
- {
- mexPrintf("%d ",j);
- mexPrintf("ija=%d ",ija_d[j]);
- mexPrintf("sa=%f\n",sa_d[j]);
- }*/
- //ija_p=(Vec_INT*)mxMalloc(sizeof(Vec_INT));
- ija_p.Format(ija_d,Size_SparseMatrixRow+1);
- //ija_p=;
- //sa_p=(Vec_DP*)mxMalloc(sizeof(Vec_DP));
- sa_p.Format(sa_d,Size_SparseMatrixRow+1);
- //sa_p=;
- /*ija_p=new Vec_INT(ija_d,Size_SparseMatrixRow+1);
- sa_p=new Vec_DP(sa_d,Size_SparseMatrixRow+1);*/
- //sprsprt(*sa_p, *ija_p);
- /*atimes(x,Err,0);
- for(i=0;i<Size*periods;i++)
- {
- Err[i]-=b[i];
- //mexPrintf("Err[%d]=%f\n",i,Err[i]);
- }*/
- BiCGStab(b,x,ITOL,TOL,ITMAX,sub_iter,err, a, indx, periods);
- mexPrintf("Estimated error: %f\n",err);
- mexPrintf("Iterations needed: %d\n",sub_iter);
- //mexPrintf("Solution vector:\n");
- for (i=0;i<y_size*(periods+y_kmin);i++)
- ya[i]=y[i];
- slowc_save=slowc;
- //mexPrintf("slowc=%f\n",slowc);
- for(i=0;i<Size;i++)
- for(j=y_kmin;j<periods+y_kmin;j++)
- {
- direction[index_vara[i]+j*y_size]=x[i+(j-y_kmin)*Size]-y[index_vara[i]+j*y_size];
- y[index_vara[i]+j*y_size]+=slowc*direction[index_vara[i]+j*y_size];
- //mexPrintf("x[%d]=y[%d]=%f\n",i+(j-y_kmin)*Size, index_vara[i]+j*y_size, y[index_vara[i]+j*y_size]);
- }
- /*if(ija_p)
- mxFree(ija_p);
- if(sa_p)
- mxFree(sa_p);*/
- time00=clock();
- if (print_it)
- {
- //pctimer_t t2 = pctimer();
- clock_t t2 = clock();
- mexPrintf("(** %f milliseconds **)\n", 1000.0*(double(t2) - double(t1))/double(CLOCKS_PER_SEC));
- mexEvalString("drawnow;");
- }
-}
-
-
-
-
-void LinBCG::asolve(Vec_I_DP &b, Vec_O_DP &x, const int itrnsp, Mat_DP &a, Vec_IO_INT &indx, const int periods)
-{
- int i;
-
- int n=b.size();
- //mexPrintf("n=%d a.nrows()=%d a.ncols()=%d\n",n, a.nrows(), a.ncols());
- //DIAG(A)
- //for(i=0;i<n;i++) x[i]=((sa_p)[i] != 0.0 ? b[i]/(sa_p)[i] : b[i]);
- //Identity
- //for(i=0;i<n;i++) x[i]=b[i];
- //LU
- /*mexPrintf("x=\n");
- x.Print();*/
-
- /*mexPrintf("b=\n");
- Vec_O_DP(b).Print();*/
- /*mexPrintf("x=\n");
- x.Print();*/
- x=b;
- if(itrnsp)
- {
- //lubksb_trp(a, indx, x);
- lubksb_blck_trp(a, indx, x, periods);
- }
- else
- {
- //lubksb(a, indx, x);
- lubksb_blck(a, indx, x, periods);
- }
- /*mexPrintf("solution\n");
- x.Print();
-
- Vec_DP z(n);
-
- atimes(x,z,itrnsp);
-
- //for(i=0;i<n;i++)
- //z[i]=z[i]-b[i];
-
- mexPrintf("check\n");
- z.Print();*/
-}
-
-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;
- //mexPrintf("sprsax\n");
- int n=x.size();
- /*mexPrintf("n=%d\n",n);
- mexPrintf("ija[0]=%d\n",ija[0]);*/
- if (ija[0] != n+1)
- mexPrintf("Error: \n sprsax: mismatched vector and matrix\n");
- /*mexPrintf("ok0a\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]];
- }
- }
- /*mexPrintf("ok0b\n");*/
-}
-
-void LinBCG::sprstx(Vec_I_DP &sa, Vec_I_INT &ija, Vec_I_DP &x, Vec_O_DP &b)
-{
- int i,j,k;
- //mexPrintf("sprstx\n");
- 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::sprsin(Vec_DP s, Vec_INT ij)
-{
- sa_p.Format(s);
- ija_p.Format(ij);
-}
-
-void LinBCG::sprsin(DP *s, int *ij, int size)
-{
- sa_p.Format(size);
- ija_p.Format(size);
- for(int i=0;i<size;i++)
- {
- sa_p[i]=s[i];
- ija_p[i]=ij[i];
- }
-}
-
-
-void LinBCG::sprsin(double *a, int NbRow, int NbCol, const DP thresh)
-{
- int i,j,k;
- int n=NbRow;
- for(i=1, j=0;i<NbRow*NbCol;i++)
- {
- //mexPrintf("%d %% %d = %d a[%d]=%f\n",i,NbRow+1, (i)%(NbRow+1),i,a[i]);
- if(fabs(a[i])>=thresh && ((i)%(NbRow+1))!=0)
- j++;
- }
- //mexPrintf("j=%d\n",j);
- ija_p.Format(j+NbRow+1);
- sa_p.Format(j+NbRow+1);
- //mexPrintf("nz=%d\n",j+NbRow+1);
- for(j=0;j<n;j++)
- sa_p[j]=a[j+j*n];
- ija_p[0]=n+1;
- k=n;
- for(i=0;i<n;i++)
- {
- for(j=0;j<n;j++)
- {
- if(fabs(a[j*n+i]) >= thresh && i != j)
- {
- sa_p[++k]=a[j*n+i];
- ija_p[k]=j;
- }
- }
- ija_p[i+1]=k+1;
- }
-}
-
-void LinBCG::sprsin(Mat_DP &a, const DP thresh)
-{
- int i,j,k=0;
-
- int n=a.nrows();
- for(i=0;i<n;i++)
- for(j=0;j<n;j++)
- if(fabs(a[i][j])>=thresh && j!=i)
- k++;
- ija_p.Format(k+n+1);
- sa_p.Format(k+n+1);
- for (j=0;j<n;j++) sa_p[j]=a[j][j];
- ija_p[0]=n+1;
- k=n;
- for (i=0;i<n;i++) {
- for (j=0;j<n;j++) {
- if (fabs(a[i][j]) >= thresh && i != j) {
- //if (++k > nmax) nrerror("sprsin: sa and ija too small");
- k++;
- sa_p[k]=a[i][j];
- ija_p[k]=j;
- }
- }
- ija_p[i+1]=k+1;
- }
-}
-
-
-void LinBCG::sprsout(Mat_DP &a)
-{
- int i,j,k,l,ke;
- double maxi=0;
- int N=ija_p[0]-1;
- //int total=ija_p[N];
- //mexPrintf("sprsout\n");
- for(i=0;i<N;i++)
- {
- l=ija_p[i];
- k=ija_p[l];
- ke=ija_p[i+1];
- for(j=0;j<N;j++)
- {
- //mexPrintf("i=%d, j=%d\n",i,j);
- if(i==j)
- //mexPrintf(fmt,sa_p[i]);
- a[i][i]=sa_p[i];
- else if(j==k)
- {
- //mexPrintf(fmt,sa_p[l]);
- a[i][j]=sa_p[l];
- if(l+1<ke)
- {
- l++;
- k=ija_p[l];
- }
- }
- else
- //mexPrintf(fmt,0.0);
- a[i][j]=0;
- }
- //mexPrintf("\n");
- }
- //mexPrintf("\n");
- //mexPrintf("end of sprsout\n");
-}
-
-
-
-void LinBCG::sprsprt()
-{
- int i,j,k,l,ke;
- double maxi=0;
- int N=ija_p[0]-1;
- int total=ija_p[N];
- char fmt[12];
- memset(fmt,'\0',12);
- string pf;
- for(i=0;i<total;i++)
- if(fabs(sa_p[i])>maxi)
- maxi=fabs(sa_p[i]);
- int size=int(log10(maxi))+1;
- if(size<7)
- {
- sprintf(fmt, "%d", size+8);
- pf= fmt;
- pf.insert(0,"% ");
- pf.append(".6f ");
- pf.copy(fmt,8);
- }
- else
- strcpy(fmt,"%e ");
- for(i=0;i<N;i++)
- {
- l=ija_p[i];
- k=ija_p[l];
- ke=ija_p[i+1];
- for(j=0;j<N;j++)
- {
- if(i==j)
- mexPrintf(fmt,sa_p[i]);
- else if(j==k)
- {
- mexPrintf(fmt,sa_p[l]);
- if(l+1<ke)
- {
- l++;
- k=ija_p[l];
- }
- }
- else
- mexPrintf(fmt,0.0);
- }
- mexPrintf("\n");
- }
- mexPrintf("\n");
-}
-
-
-void LinBCG::sprs_sprt()
-{
- int i,j,k,l,ke;
- double maxi=0;
- int N=ija_p[0]-1;
- mexPrintf("N=%d\n",N);
- int total=ija_p[N];
- mexPrintf("total=%d\n",total);
- char fmt[12], fmts[12], fmtd[12];
- memset(fmt,'\0',12);
- memset(fmts,'\0',12);
- memset(fmtd,'\0',12);
- string pf, pf1, pf2;
- for(i=0;i<total;i++)
- if(fabs(sa_p[i])>maxi)
- maxi=fabs(sa_p[i]);
- int size=int(log10(maxi))+1;
- if(size<7)
- {
- sprintf(fmt, "%d", size+8);
- pf= fmt;
- pf.insert(0,"% ");
- pf1=pf2=pf;
- pf.append(".6f ");
- pf1.append("s ");
- pf2.append("d ");
- pf.copy(fmt,8);
- pf1.copy(fmts,8);
- pf2.copy(fmtd,8);
- }
- else
- strcpy(fmt,"%e ");
- mexPrintf(fmts,"i");
- mexPrintf(fmts,"ija");
- mexPrintf(fmts,"sa");
- if(ijat_p.size())
- {
- mexPrintf(fmts,"ijat");
- mexPrintf(fmts,"sat");
- }
- mexPrintf("\n");
- for(i=0;i<total;i++)
- {
- mexPrintf(fmtd,i);
- mexPrintf(fmtd,ija_p[i]);
- mexPrintf(fmt,sa_p[i]);
- if(ijat_p.size())
- {
- mexPrintf(fmtd,ijat_p[i]);
- mexPrintf(fmt,sat_p[i]);
- }
- mexPrintf("\n");
- }
- mexPrintf("\n");
-}
-
-void LinBCG::BiCG(Vec_I_DP &b, Vec_IO_DP &x, const int itol, const DP tol,
- const int itmax, int &iter, DP &err, Mat_DP &a, Vec_IO_INT &indx, const int periods)
-{
- DP ak,akden,bk,bkden=1.0,bknum,bnrm=0,dxnrm,xnrm,zm1nrm,znrm=0;
- const DP EPS=1.0e-14;
- int j;
- int n=b.size();
- //mexPrintf("n=%d\n",n);
- 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, a, indx, periods);
- }
- else if (itol == 2)
- {
- asolve(b,z,0, a, indx, periods);
- bnrm=snrm(z,itol);
- asolve(r,z,0, a, indx, periods);
- }
- else if (itol == 3 || itol == 4)
- {
- asolve(b,z,0, a, indx, periods);
- bnrm=snrm(z,itol);
- asolve(r,z,0, a, indx, periods);
- znrm=snrm(z,itol);
- }
- else
- mexPrintf("Error: \n illegal itol in BiCG\n");
- //cout << fixed << setprecision(6);
- while (iter < itmax)
- {
- ++iter;
- asolve(rr,zz,1, a, indx, periods);
- 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, a, indx, periods);
- 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;
- mexPrintf("iter=%d err=%f ak=%f akden=%f\n",iter,err,ak,akden);
- if (err <= tol)
- break;
- }
-}
-
-void LinBCG::BiCGStab(Vec_I_DP &b, Vec_IO_DP &x, const int itol, const DP tol,
- const int itmax, int &iter, DP &err, Mat_DP &a, Vec_IO_INT &indx, const int periods)
-{
- DP akden,bkden=1.0,bknum,bnrm=0,dxnrm,xnrm,zm1nrm,znrm=0, alpha, omega, beta, omega_num, omega_den, snorm;
- const DP EPS=1.0e-14;
- int j;
- int n=b.size();
- //mexPrintf("n=%d\n",n);
- Vec_DP p(n),v(n),r(n),rr(n),z(n),s(n), t(n), phat(n), shat(n);
- iter=0;
- // Compute the residual
- 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);
- //mexPrintf("asolve(r,z,..)\n");
- asolve(r,z,0, a, indx, periods);
- }
- else if (itol == 2)
- {
- asolve(b,z,0, a, indx, periods);
- bnrm=snrm(z,itol);
- asolve(r,z,0, a, indx, periods);
- }
- else if (itol == 3 || itol == 4)
- {
- asolve(b,z,0, a, indx, periods);
- bnrm=snrm(z,itol);
- asolve(r,z,0, a, indx, periods);
- znrm=snrm(z,itol);
- }
- else
- mexPrintf("Error: \n illegal itol in BiCGStab\n");
- //cout << fixed << setprecision(6);
- while (iter < itmax)
- {
- ++iter;
- //asolve(rr,zz,1, a, indx);
- for (bknum=0.0,j=0;j<n;j++)
- bknum += r[j]*rr[j];
- if (iter == 1)
- {
- for (j=0;j<n;j++)
- {
- p[j]=r[j];
- }
- }
- else
- {
- beta=bknum/bkden*(alpha/omega);
- for (j=0;j<n;j++)
- {
- p[j]=r[j]+beta*(p[j]-omega*v[j]);
- }
- }
- //mexPrintf("asolve(p,phat,..)\n");
- asolve(p,phat,1, a, indx, periods);
-
- atimes(phat,v,0);
- for(akden=0.0,j=0;j<n;j++)
- akden += rr[j]*v[j];
- alpha=bknum/akden;
- for(snorm=0.0, j=0;j<n;j++)
- {
- s[j]=r[j]-alpha*v[j];
- snorm+=s[j]*s[j];
- }
- if(sqrt(snorm)<tol)
- {
- for(j=0;j<n;j++)
- x[j]+=alpha*phat[j]+omega*shat[j];
- break;
- }
- //mexPrintf("asolve(s,shat,..)\n");
- asolve(s,shat,1, a, indx, periods);
- atimes(shat,t,0);
- for(omega_num=0.0,omega_den=0.0,j=0;j<n;j++)
- {
- omega_num+=t[j]*s[j];
- omega_den+=t[j]*t[j];
- }
- omega=omega_num/omega_den;
- for(j=0;j<n;j++)
- {
- x[j]+=alpha*phat[j]+omega*shat[j];
- r[j]=s[j]-omega*t[j];
- }
- bkden=bknum;
-
-
- 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(alpha)*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;
- mexPrintf("iter=%d err=%f alpha=%f akden=%f\n",iter,err,alpha,akden);
- mexEvalString("drawnow;");
- if (err <= tol)
- break;
- }
-}
-
-
-void LinBCG::sprs_col_index()
-{
- int i, j, k, begin, end;
- int nze=sa_p.size();
- ijat_p.Format(nze);
- sat_p.Format(nze);
- begin=ija_p[0];
- int n=begin-1;
- //mexPrintf("n=%d\n",n);
- map<std::pair<int, int>, double> list;
- map<std::pair<int, int>, double>::iterator list_it;
- for(j=0;j<n;j++)
- sat_p[j]=sa_p[j];
- sat_p[n]=sa_p[n];
- ijat_p[n]=ija_p[n];
- for(i=0;i<n;i++)
- {
- end=ija_p[i+1];
- for(k=begin;k<end;k++)
- {
- //mexPrintf("make_pair(ija_p[k]=%d, i=%d)=%f\n",ija_p[k],i,sa_p[k]);
- list[make_pair(ija_p[k],i)]=sa_p[k];
- }
- begin=end;
- }
- int col=-1;
- int cur_pos=n+1;
- list_it=list.begin();
- while(list_it!=list.end())
- {
- //mexPrintf("list_it->first.first=%d, list_it->first.second=%d, list_it->second=%f\n",list_it->first.first, list_it->first.second, list_it->second);
- while(col!=list_it->first.first)
- {
- //while(col!=list_it->first.first
- col++;
- ijat_p[col]=cur_pos;
- }
- ijat_p[cur_pos]=list_it->first.second;
- sat_p[cur_pos++]=list_it->second;
- list_it++;
- }
-}
-
-
-void LinBCG::sprs_swap_line_copy(map<std::pair<int, int>, double> &list, int &pos_2_blck, int begin, int end)
-{
- int i, j, init_pos;
- for(i=begin;i<end;i++)
- {
- init_pos=pos_2_blck;
- list[make_pair(i,init_pos)]=sa_p[i];
- for(j=ija_p[i];j<ija_p[i+1];j++)
- {
- //mexPrintf("list[%d,%d]=%f\n",pos_2_blck,ija_p[j],sa_p[j]);
- list[make_pair(pos_2_blck++,ija_p[j])]=sa_p[j];
- }
- }
-}
-
-
-void LinBCG::sprs_swap_line_exchange(map<std::pair<int, int>, double> &list, int &pos_2_blck, int LS, int LD)
-{
- int j;
- bool OK=false, OK1=true;
- int init_pos=pos_2_blck;
- for(j=ija_p[LS];j<ija_p[LS+1];j++)
- {
- if(ija_p[j]==LD)
- {
- OK=true;
- list[make_pair(LD,init_pos)]=sa_p[j];
- //mexPrintf("found list[%d,%d]=%f\n",LD,init_pos,sa_p[j]);
- }
- else
- {
- if(ija_p[j]>=LS && OK1)
- {
- OK1=false;
- if(sa_p[LS]!=0.0)
- {
- //mexPrintf("list[%d,%d]=%f\n",pos_2_blck,LS,sa_p[LS]);
- list[make_pair(pos_2_blck++,LS)]=sa_p[LS];
- }
- }
- //mexPrintf("list[%d,%d]=%f\n",pos_2_blck,ija_p[j],sa_p[j]);
- list[make_pair(pos_2_blck++,ija_p[j])]=sa_p[j];
- }
- }
- if(OK1)
- {
- if(sa_p[LS]!=0.0)
- {
- //mexPrintf("list[%d,%d]=%f\n",pos_2_blck,LS,sa_p[LS]);
- list[make_pair(pos_2_blck++,LS)]=sa_p[LS];
- }
- }
- if(!OK)
- {
- list[make_pair(LD,init_pos)]=0.0;
- //mexPrintf("not found list[%d,%d]=%f\n",LD,init_pos,0.0);
- }
-}
-
-void LinBCG::sprs_swap_line(int L0, int L1)
-{
- int i,j,k, pos_2_blck, init_pos;
- map<std::pair<int, int>, double> list;
- map<std::pair<int, int>, double>::iterator list_it;
- int n=ija_p[0]-1;
- if(L0>L1)
- {
- i=L0;
- L0=L1;
- L1=i;
- }
- else if(L0==L1)
- return;
- else if(L0>=n || L1>=n)
- {
- mexPrintf("out of range in sprs_swap_line (L0=%d and L1=%d)\n",L0, L1);
- mexEvalString("drawnow;");
- mexEvalString("st=fclose('all');clear all;");
- filename+=" stopped";
- mexErrMsgTxt(filename.c_str());
- }
- int nze=sa_p.size();
- bool OK, OK1;
- int nze0=0, nze1=0;
- pos_2_blck=n+1;
-
- sprs_swap_line_copy(list, pos_2_blck, 0, L0);
-
- //mexPrintf("L1=>L0\n");
-
- sprs_swap_line_exchange(list, pos_2_blck, L1, L0);
-
- sprs_swap_line_copy(list, pos_2_blck, L0+1, L1);
-
- //mexPrintf("L0=>L1\n");
-
- sprs_swap_line_exchange(list, pos_2_blck, L0, L1);
-
- sprs_swap_line_copy(list, pos_2_blck, L1+1, n);
-
- sa_p.Format(list.size()+1);
- ija_p.Format(list.size()+1);
- //mexPrintf("list.size()+1=%d\n",list.size()+1);
- list_it=list.begin();
- while(list_it!=list.end())
- {
- //mexPrintf("sa_p[list_it->first.first=%d]=list_it->second=%f\n",list_it->first.first,list_it->second);
- sa_p[list_it->first.first]=list_it->second;
- //mexPrintf("ija_p[list_it->first.first=%d]=list_it->first.second=%d\n",list_it->first.first,list_it->first.second);
- ija_p[list_it->first.first]=list_it->first.second;
- list_it++;
- }
- ija_p[n]=list.size()+1;
- sa_p[n]=0;
- sprs_col_index();
-}
-
-
-
-void LinBCG::sprsludcmp(Mat_DP &a, Vec_O_INT &indx, DP &d)
-{
- const DP TINY=1.0e-20;
- int i,imax,j,k, ii, ik, begin_ija, end_ija, begin_ijat, end_ijat, dum1, dum2;
- DP big,dum,sum,temp;
- Vec_DP saa_p(sa_p);
- Vec_INT ijaa_p(ija_p);
- int nze=sa_p.size();
- int n=ija_p[0]-1;
- Vec_DP vv(n);
- d=1.0;
-
- n=a.nrows();
- //a=(double**)mxMalloc(n*n*sizeof(double));
- //mexPrintf("sprsludcmp\n");
-
- for (i=0;i<n;i++) {
- big=0.0;
- for (j=0;j<n;j++)
- {
- //mexPrintf("%f ",a[i][j]);
- if ((temp=fabs(a[i][j])) > big) big=temp;
- }
- //mexPrintf("\n");
- if(big == 0.0)
- {
- mexPrintf("Singular Matrix in routine sprsludcmp\n");
- mexEvalString("st=fclose('all');clear all;");
- mexErrMsgTxt("End of simulate");
- }
- vv[i]=1.0/big;
- }
- /*mexPrintf("First Step\n");
- mexPrintf("a.ncols()=%d\n",a.ncols());
- mexPrintf("a.nrows()=%d\n",a.nrows());*/
- /*for(i=0;i<n;i++)
- {
- big=fabs(sa_p[i]);
- for(j=ija_p[i];j<ija_p[i+1];j++)
- if(temp=fabs(sa_p[j])>big)
- big=temp;
- if(big==0.0)
- {
- mexPrintf("Singular Matrix in routine sprsludcmp\n");
- mexErrMsgTxt("End of simulate");
- }
- vv[i]=1.0/big;
- }*/
- /*sprs_col_index(); // create a col-index compact storage
- for (j=0;j<n;j++)
- {
- //imax=j;
- for (i=ijat_p[j];i<ijat_p[j+1] && ijat_p[i]<j;i++)
- {
- sum=sat_p[i];
- begin_ija=ija_p[i];
- end_ija=ija_p[i+1];
- begin_ijat=ijat_p[i];
- end_ijat=ijat_p[i+1];
- while(begin_ija<end_ija && begin_ijat<end_ijat && begin_ijat<begin_ija)
- {
- if(ija_p[begin_ija]<ijat_p[begin_ijat])
- begin_ija++;
- else if(ija_p[begin_ija]>ijat_p[begin_ijat])
- begin_ijat++;
- else
- sum -= sa_p[begin_ija++]*sat_p[begin_ijat++];
- }
- sat_p[i]=sum;
- }
- big=0.0;
- for (;i<ijat_p[j+1] ;i++)
- {
- sum=sat_p[i];
- begin_ija=ija_p[i];
- end_ija=ija_p[i+1];
- begin_ijat=ijat_p[i];
- end_ijat=ijat_p[i+1];
- while(begin_ija<end_ija && ija_p[begin_ija]<j && begin_ijat<end_ijat && ijat_p[begin_ijat]<j)
- {
- if(ija_p[begin_ija]<ijat_p[begin_ijat])
- begin_ija++;
- else if(ija_p[begin_ija]>ijat_p[begin_ijat])
- begin_ijat++;
- else
- sum -= sa_p[begin_ija++]*sat_p[begin_ijat++];
- }
- sat_p[i]=sum;
- if((dum=vv[ija_p[i]]*fabs(sum)) >= big)
- {
- big=dum;
- mexPrintf("imax=i=%d\n",i);
- imax=ija_p[i];
- }
- }
- if (j != imax)
- {
- sprs_swap_line(j,imax);
- //for (k=0;k<n;k++)
- // {
- // dum=a[imax][k];
- // a[imax][k]=a[j][k];
- // a[j][k]=dum;
- // }
- d = -d;
- vv[imax]=vv[j];
- }
- indx[j]=imax;
- if (sa_p[j] == 0.0) sa_p[j]=TINY;
- if (j != n-1)
- {
- dum=1.0/(sat_p[j]);
- for (i=ija_p[j];i<ija_p[j+1];i++)
- if(ija_p[j]>j)
- sa_p[i] *= dum;
- }
- }
-}*/
- for (j=0;j<n;j++)
- {
- for (i=0;i<j;i++)
- {
- sum=a[i][j];
- for (k=0;k<i;k++)
- sum -= a[i][k]*a[k][j];
- a[i][j]=sum;
- }
- big=0.0;
- //mexPrintf("ok0\n");
- for (i=j;i<n;i++)
- {
- sum=a[i][j];
- for (k=0;k<j;k++)
- sum -= a[i][k]*a[k][j];
- a[i][j]=sum;
- if ((dum=vv[i]*fabs(sum)) >= big)
- {
- big=dum;
- imax=i;
- }
- }
- //mexPrintf("ok1\n");
- if (j != imax)
- {
- for (k=0;k<n;k++)
- {
- dum=a[imax][k];
- a[imax][k]=a[j][k];
- a[j][k]=dum;
- }
- d = -d;
- vv[imax]=vv[j];
- }
- indx[j]=imax;
- //mexPrintf("ok2\n");
- if (a[j][j] == 0.0)
- a[j][j]=TINY;
- if (j != n-1)
- {
- dum=1.0/(a[j][j]);
- for (i=j+1;i<n;i++)
- a[i][j] *= dum;
- }
- //mexPrintf("ok3\n");
- }
- mexPrintf("end of sprsludcmp\n");
- //sprsin(a,1e-10);
- //mxFree(a);
-}
-
-
-
-void LinBCG::lubksb_blck(Mat_I_DP &a, Vec_I_INT &indx, Vec_IO_DP &b, const int periods)
-{
- int i,ii,ip,j, k;
- DP sum;
- int n=a.nrows();
- for(k=0;k<periods;k++)
- {
- ii=0;
- for (i=0;i<n;i++)
- {
- ip=indx[i];
- sum=b[ip+k*n];
- b[ip+k*n]=b[i+k*n];
- /*if (ii != 0)
- for (j=ii-1;j<i;j++)
- sum -= a[i][j]*b[j+k*n];
- else if (sum != 0.0)
- ii=i+1;*/
- for (j=0;j<i;j++)
- sum -= a[i][j]*b[j+k*n];
- b[i+k*n]=sum;
- }
- }
- for(k=periods-1;k>=0;k--)
- {
- for (i=n-1;i>=0;i--)
- {
- sum=b[i+k*n];
- for (j=i+1;j<n;j++)
- sum -= a[i][j]*b[j+k*n];
- if(fabs(a[i][i])<tol)
- {
- mexPrintf("Error: Singularity in lubksb_blck a[%d][%d](max=%d)=%g\n",i,i,n,a[i][i]);
- mexEvalString("st=fclose('all');clear all;");
- mexErrMsgTxt("Exit from Dynare");
- }
- b[i+k*n]=sum/a[i][i];
- }
- }
-}
-
-void LinBCG::lubksb_blck_trp(Mat_I_DP &a, Vec_I_INT &indx, Vec_IO_DP &b, const int periods)
-{
- int i,ii,ip,j, k;
- DP sum;
- int n=a.nrows();
- for(k=0;k<periods;k++)
- {
- ii=0;
- for (i=0;i<n;i++)
- {
- ip=indx[i];
- sum=b[ip+k*n];
- b[ip+k*n]=b[i+k*n];
- /*if (ii != 0)
- for (j=ii-1;j<i;j++)
- sum -= a[j][i]*b[j+k*n];
- else if (sum != 0.0)
- ii=i+1;*/
- for (j=0;j<i;j++)
- sum -= a[j][i]*b[j+k*n];
- if(fabs(a[i][i])<tol)
- {
- mexPrintf("Error: Singularity in lubksb_blck_trp a[%d][%d](max=%d)=%g\n",i,i,n,a[i][i]);
- mexEvalString("st=fclose('all');clear all;");
- mexErrMsgTxt("Exit from Dynare");
- }
- b[i+k*n]=sum/a[i][i];
- }
- }
- for(k=periods-1;k>=0;k--)
- {
- for (i=n-1;i>=0;i--)
- {
- sum=b[i+k*n];
- for (j=i+1;j<n;j++)
- sum -= a[j][i]*b[j+k*n];
- b[i+k*n]=sum;
- }
- }
-}
-
-
-
-void LinBCG::lubksb(Mat_I_DP &a, Vec_I_INT &indx, Vec_IO_DP &b)
-{
- // Solving A.x=b
- int i,ii,ip,j;
- DP sum;
- int n=a.nrows();
- ii=0;
- // Solving L.y=b
- // L is a lower triangular matrix with element on the main diagonal equal to 1
- // Implementation of a foreward substitution to solve this system
- for (i=0;i<n;i++)
- {
- ip=indx[i];
- sum=b[ip];
- b[ip]=b[i];
- if (ii != 0)
- for (j=ii-1;j<i;j++)
- sum -= a[i][j]*b[j];
- else if (sum != 0.0)
- ii=i+1;
- b[i]=sum;
- }
- // Solving U.x=y
- // U is an upper triangular matrix with element on the main diagonal is different from 1
- // Implementation of a backward substitution to solve this system
- for (i=n-1;i>=0;i--)
- {
- sum=b[i];
- for (j=i+1;j<n;j++)
- sum -= a[i][j]*b[j];
- if(fabs(a[i][i])<tol)
- {
- mexPrintf("Error: Singularity in lubksb a[%d][%d](max=%d)=%g\n",i,i,n,a[i][i]);
- mexEvalString("st=fclose('all');clear all;");
- mexErrMsgTxt("Exit from Dynare");
- }
- b[i]=sum/a[i][i];
- }
-}
-
-
-void LinBCG::lubksb_trp(Mat_I_DP &a, Vec_I_INT &indx, Vec_IO_DP &b)
-{
- // Solving A'.x=b <=> (LU)'.x=b <=> U'L'x=b <=> U'.y=b with x solution of L'x=y
- int i,ii,ip,j;
- DP sum;
- int n=a.nrows();
- ii=0;
- // Solving U'.y=b
- // U' is a lower triangular matrix with element on the main diagonal differnet from 1
- // Implementation of a foreward substitution to solve this system
- for (i=0;i<n;i++)
- {
- ip=indx[i];
- sum=b[ip];
- b[ip]=b[i];
- if (ii != 0)
- for (j=ii-1;j<i;j++)
- sum -= a[j][i]*b[j];
- else if (sum != 0.0)
- ii=i+1;
- if(fabs(a[i][i])<tol)
- {
- mexPrintf("Error: Singularity in lubksb_trp a[%d][%d](max=%d)=%g\n",i,i,n,a[i][i]);
- mexEvalString("st=fclose('all');clear all;");
- mexErrMsgTxt("Exit from Dynare");
- }
- b[i]=sum/a[i][i];
- }
- // Solving L'.x=y
- // L' is an upper triangular matrix with element on the main diagonal equal to 1
- // Implementation of a backward substitution to solve this system
- for (i=n-1;i>=0;i--)
- {
- sum=b[i];
- for (j=i+1;j<n;j++)
- sum -= a[j][i]*b[j];
- b[i]=sum;
- }
-}
-
-
-void LinBCG::sprs_LU(Vec_O_DP &sL, Vec_O_INT &ijL, Vec_O_DP &sU, Vec_O_INT &ijU)
-{
- int j,k, nU, nL;
- int n=ija_p[0]-1;
- int L_count, A_count, U_count;
- L_count=A_count=U_count=n+1;
- int nn=ija_p.size();
- Vec_INT indx(n);
- Vec_INT save_ija_p(ija_p);
- Vec_DP save_sa_p(sa_p);
- DP d;
-
- mexPrintf("sprs_LU\n");
-
- Mat_DP a(n,n);
-
- sprsludcmp(a,indx, d);
-
- Vec_DP res(n), res_save(n);
- for(int i=0;i<n;i++)
- res_save[i]=res[i]=i;
-
- a.Print();
-
- //lubksb(a, indx, res, 0);
- lubksb(a, indx, res);
-
- mexPrintf("solution\n");
- res.Print();
-
- Vec_DP z(n);
-
- atimes(res,z,0);
-
- mexPrintf("check\n");
- z.Print();
-
- nU=nL=n;
- for(k=0;k<n;k++)
- for(j=ija_p[k];j<ija_p[k+1];j++)
- if(ija_p[j]<k)
- nU++;
- else
- nL++;
- mexPrintf("nU+1=%d, nL+1=%d sa_p.size()=%d\n",nU+1, nL+1, sa_p.size());
- ijU.Format(nU+1);
- sU.Format(nU+1);
-
- ijL.Format(nL+1);
- //mexPrintf("ijL[0]=%d\n",ijL[0]);
- sL.Format(nL+1);
-
- for(k=0;k<n;k++)
- {
- //the diagonal elements
- sU[k]=sa_p[k];
- sL[k]=1.0;
-
- ijU[k]=U_count;
- ijL[k]=L_count;
- for(j=ija_p[k];j<ija_p[k+1];j++)
- {
- if(ija_p[j]<k)
- {
- mexPrintf("U_count=%d, A_count=%d\n",U_count, A_count);
- ijU[U_count]=ija_p[A_count];
- sU[U_count++]=sa_p[A_count++];
- }
- else
- {
- mexPrintf("L_count=%d, A_count=%d\n",L_count, A_count);
- ijL[L_count]=ija_p[A_count];
- sL[L_count++]=sa_p[A_count++];
- }
- }
- }
- sU[n]=0.0;
- sL[n]=0.0;
- ijU[n]=nU+1;
- ijL[n]=nL+1;
- mexPrintf("ijL[0]=%d\n",ijL[0]);
- mexPrintf("end0 of sprs_LU\n");
- sa_p.Format(save_sa_p);
- ija_p.Format(save_ija_p);
- mexPrintf("end of sprs_LU\n");
-}
-
-void LinBCG::Initialize(string filename_arg, double res1_arg, double res2_arg, double max_res_arg, double slowc_arg, double *ya_arg, double *direction_arg, int iter_arg)
-{
- filename=filename_arg;
- res1=res1_arg;
- res2=res2_arg;
- max_res=max_res_arg;
- slowc=slowc_arg;
- slowc_save=slowc;
- ya=ya_arg;
- direction=direction_arg;
- iter=iter_arg;
-}
diff --git a/mex/sources/simulate/linbcg.hh b/mex/sources/simulate/linbcg.hh
deleted file mode 100644
index b040c403d22a13f877633bffa845d10b668a71e7..0000000000000000000000000000000000000000
--- a/mex/sources/simulate/linbcg.hh
+++ /dev/null
@@ -1,458 +0,0 @@
-/*
- * Copyright (C) 2007-2008 Dynare Team
- *
- * This file is part of Dynare.
- *
- * Dynare is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * Dynare is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with Dynare. If not, see <http://www.gnu.org/licenses/>.
- */
-
-#ifndef LINBCG_HH_INCLUDED
-#define LINBCG_HH_INCLUDED
-
-#include <iostream>
-#include <iomanip>
-#include <cmath>
-#include <complex>
-#include <map>
-#include <string>
-#include <ctime>
-#include "mex.h"
-
-using namespace std;
-
-const int debile=10;
-const double tol=1e-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
- // NRVec & operator=(const T* rhs);
- inline T & operator[](const int i); //i'th element
- inline const T & operator[](const int i) const;
- inline int size() const;
- inline void Format(T* rhs, int n);
- inline void Format(const NRVec rhs);
- inline void Format(int n);
- inline void Print();
- ~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]*/(T*)mxMalloc(sizeof(T)*n)) {}
-
-template <class T>
-NRVec<T>::NRVec(const T& a, int n) : nn(n), v(/*new T[n]*/(T*)mxMalloc(sizeof(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]*/(T*)mxMalloc(sizeof(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]*/(T*)mxMalloc(sizeof(T)*rhs.nn))
-{
- for(int i=0; i<nn; i++)
- v[i] = rhs[i];
-}
-
-template <class T>
-void NRVec<T>::Format(T* rhs, int n)
-{
- if(v)
- mxFree(v);
- v=NULL;
- nn=n;
- v=(T*)mxMalloc(sizeof(T)*n);
- for(int i=0; i<nn; i++)
- v[i] = rhs[i];
-}
-
-template <class T>
-void NRVec<T>::Format(const NRVec rhs)
-{
- if(v)
- mxFree(v);
- v=NULL;
- nn=rhs.nn;
- v=(T*)mxMalloc(sizeof(T)*nn);
- for(int i=0; i<nn; i++)
- v[i] = rhs[i];
-}
-
-
-template <class T>
-void NRVec<T>::Format(int n)
-{
- if(v)
- mxFree(v);
- v=NULL;
- nn=n;
- v=(T*)mxMalloc(sizeof(T)*n);
-}
-
-
-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)
- mxFree(v);
- nn=rhs.nn;
- v= (T*)mxMalloc(sizeof(T)*nn);
- }
- for (int i=0; i<nn; i++)
- v[i]=rhs[i];
- }
- return *this;
-}
-
-template <class T>
-void NRVec<T>::Print()
-{
- int i;
- for(i=0;i<nn;i++)
- mexPrintf("%f\n",v[i]);
-}
-
-/*template <class T>
-NRVec<T> & NRVec<T>::operator=(const 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 (v != 0)
- mxFree(v);
- v= (T*)mxMalloc(sizeof(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
-{
- if(i>=nn || i<0)
- mexPrintf("Error: out of range in Vect[] operator\n");
- return v[i];
-}
-
-template <class T>
-inline const T & NRVec<T>::operator[](const int i) const //subscripting
-{
- if(i>=nn || i<0)
- mexPrintf("Error: out of range in Vect[] operator\n");
- 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)
- mxFree(v);
-}
-
-
-template <class T>
-class NRMat {
-private:
- int nn;
- int mm;
- T **v;
-public:
- NRMat();
- NRMat(int n, int m); // Zero-based array
- NRMat(const T &a, int n, int m); //Initialize to constant
- NRMat(const T *a, int n, int m); // Initialize to array
- NRMat(const NRMat &rhs); // Copy constructor
- NRMat & operator=(const NRMat &rhs); //assignment
- NRMat & operator=(const T &a); //assign a to every element
- inline T* operator[](const int i); //subscripting: pointer to row i
- inline const T* operator[](const int i) const;
- inline int nrows() const;
- inline int ncols() const;
- inline void Format(const int n,const int m);
- void Print();
- ~NRMat();
-};
-
-template <class T>
-NRMat<T>::NRMat() : nn(0), mm(0), v(0) {}
-
-template <class T>
-NRMat<T>::NRMat(int n, int m) :
- nn(n),
- mm(m),
- v((T**)mxMalloc(sizeof(T*)*n))
-{
- v[0] = /*new T[m*n]*/(T*)mxMalloc(sizeof(T)*m*n);
- for (int i=1; i< n; i++)
- v[i] = v[i-1] + m;
-}
-
-template <class T>
-NRMat<T>::NRMat(const T &a, int n, int m) : nn(n), mm(m), v(/*new T*[n]*/(T**)mxMalloc(sizeof(T*)*n))
-{
- int i,j;
- v[0] = /*new T[m*n]*/(T*)mxMalloc(sizeof(T)*m*n);
- for (i=1; i< n; i++)
- v[i] = v[i-1] + m;
- for (i=0; i< n; i++)
- for (j=0; j<m; j++)
- v[i][j] = a;
-}
-
-template <class T>
-NRMat<T>::NRMat(const T *a, int n, int m) : nn(n), mm(m), v(/*new T*[n]*/(T**)mxMalloc(sizeof(T*)*n))
-{
- int i,j;
- v[0] = /*new T[m*n]*/(T*)mxMalloc(sizeof(T)*m*n);
- for (i=1; i< n; i++)
- v[i] = v[i-1] + m;
- for (i=0; i< n; i++)
- for (j=0; j<m; j++)
- v[i][j] = *a++;
-}
-
-template <class T>
-NRMat<T>::NRMat(const NRMat &rhs) : nn(rhs.nn), mm(rhs.mm), v(/*new T*[nn]*/(T**)mxMalloc(sizeof(T*)*nn))
-{
- int i,j;
- v[0] = /*new T[mm*nn]*/(T*)mxMalloc(sizeof(T)*mm*nn);
- for (i=1; i< nn; i++)
- v[i] = v[i-1] + mm;
- for (i=0; i< nn; i++)
- for (j=0; j<mm; j++)
- v[i][j] = rhs[i][j];
-}
-
-template <class T>
-NRMat<T> & NRMat<T>::operator=(const NRMat<T> &rhs)
-// postcondition: normal assignment via copying has been performed;
-// if matrix and rhs were different sizes, matrix
-// has been resized to match the size of rhs
-{
- if (this != &rhs) {
- int i,j;
- if (nn != rhs.nn || mm != rhs.mm) {
- if (v != 0)
- {
- mxFree(v[0]);
- mxFree(v);
- /*delete[] (v[0]);
- delete[] (v);*/
- }
- nn=rhs.nn;
- mm=rhs.mm;
- //v = new T*[nn];
- v = (T**)mxMalloc(sizeof(T*)*nn);
- v[0] = /*new T[mm*nn]*/(T*)mxMalloc(sizeof(T)*mm*nn);
- }
- for (i=1; i< nn; i++)
- v[i] = v[i-1] + mm;
- for (i=0; i< nn; i++)
- for (j=0; j<mm; j++)
- v[i][j] = rhs[i][j];
- }
- return *this;
-}
-
-template <class T>
-NRMat<T> & NRMat<T>::operator=(const T &a) //assign a to every element
-{
- for (int i=0; i< nn; i++)
- for (int j=0; j<mm; j++)
- v[i][j] = a;
- return *this;
-}
-
-template <class T>
-inline T* NRMat<T>::operator[](const int i) //subscripting: pointer to row i
-{
- if(i>=nn || i<0)
- mexPrintf("Error: out of range in Mat[] operator\n");
- return v[i];
-}
-
-template <class T>
-inline const T* NRMat<T>::operator[](const int i) const
-{
- if(i>=nn || i<0)
- mexPrintf("Error: out of range in Mat[] operator\n");
- return v[i];
-}
-
-template <class T>
-inline int NRMat<T>::nrows() const
-{
- return nn;
-}
-
-template <class T>
-inline int NRMat<T>::ncols() const
-{
- return mm;
-}
-
-template <class T>
-inline void NRMat<T>::Format(const int n,const int m)
-{
- if (v != 0)
- {
- /*delete[] (v[0]);
- delete[] (v);*/
- mxFree(v[0]);
- mxFree(v);
- }
- nn=n;
- mm=m;
- v = (T**)mxMalloc(sizeof(T*)*n);//new T*[n];
- v[0] = (T*)mxMalloc(sizeof(T)*m*n);//new T[m*n];
- for (int i=1; i< n; i++)
- v[i] = v[i-1] + m;
-}
-
-
-
-template <class T>
-NRMat<T>::~NRMat()
-{
- if (v != 0)
- {
- mxFree(v[0]);
- mxFree(v);
- /*delete[] (v[0]);
- delete[] (v);*/
- }
-}
-
-
-template <class T>
-void NRMat<T>::Print()
-{
- int i,j;
- mexPrintf("nn=%d, mm=%d\n",nn,mm);
- for(i=0;i<nn;i++)
- {
- for(j=0;j<mm;j++)
- mexPrintf("%f ",v[i][j]);
- mexPrintf("\n");
- }
-}
-
-
-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;
-
-typedef const NRMat<DP> Mat_I_DP;
-typedef NRMat<DP> Mat_DP, Mat_O_DP, Mat_IO_DP;
-
-
-
-
-class LinBCG
-{
- Vec_INT ija_p, ijat_p;
- Vec_DP sa_p, sat_p;
- private:
-
-
- clock_t time00;
- string filename;
- double res1, res2, slowc, slowc_save, *ya, *direction, max_res;
- int iter;
-
- void sprs_swap_line_copy(map<std::pair<int, int>, double> &list, int &pos_2_blck, int begin, int end);
- void sprs_swap_line_exchange(map<std::pair<int, int>, double> &list, int &pos_2_blck, int LS, int LD);
-
- public:
- LinBCG();
- void SolveLinear(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, int y_size, double* y, bool print_it, bool cvg, Mat_DP &a, Vec_IO_INT &indx);
- void Preconditioner(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, int y_size, double* y, bool print_it, int type, Mat_O_DP &a, Vec_O_INT &indx);
- void asolve(Vec_I_DP &b, Vec_O_DP &x, const int itrnsp, Mat_DP &a, Vec_IO_INT &indx, const int periods);
- 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 sprsin(double *a, int NbRow, int NbCol, const DP thresh);
- void sprsin(Vec_DP s, Vec_INT ij);
- void sprsin(DP *s, int *ij, int size);
- void sprsin(Mat_DP &a, const DP thresh);
- void sprsout(Mat_DP &a);
- void BiCG(Vec_I_DP &b, Vec_IO_DP &x, const int itol, const DP tol,
- const int itmax, int &iter, DP &err, Mat_DP &a, Vec_IO_INT &indx, const int periods);
- void BiCGStab(Vec_I_DP &b, Vec_IO_DP &x, const int itol, const DP tol,
- const int itmax, int &iter, DP &err, Mat_DP &a, Vec_IO_INT &indx, const int periods);
- void sprsprt();
- void sprs_sprt();
- void sprs_col_index();
- void sprs_swap_line(int L0, int L1);
- void sprsludcmp(Mat_DP &a, Vec_O_INT &indx, DP &d);
- void lubksb_blck(Mat_I_DP &a, Vec_I_INT &indx, Vec_IO_DP &b, const int periods);
- void lubksb_blck_trp(Mat_I_DP &a, Vec_I_INT &indx, Vec_IO_DP &b, const int periods);
- void lubksb(Mat_I_DP &a, Vec_I_INT &indx, Vec_IO_DP &b);
- void lubksb_trp(Mat_I_DP &a, Vec_I_INT &indx, Vec_IO_DP &b);
- void sprs_LU(Vec_O_DP &sL, Vec_O_INT &ijL, Vec_O_DP &sU, Vec_O_INT &ijU);
- void Initialize(string filename_arg, double res1_arg, double res2_arg, double max_res_arg, double slowc_arg, double *ya_arg, double *direction_arg, int iter_arg);
-};
-
-
-#endif
diff --git a/mex/sources/simulate/simulate.cc b/mex/sources/simulate/simulate.cc
index 931f16ef3e9237485c8009a2cfa53245353efbe7..d176900f64296033417ac32b7c8966e7c2b5e53c 100644
--- a/mex/sources/simulate/simulate.cc
+++ b/mex/sources/simulate/simulate.cc
@@ -28,7 +28,6 @@
#include "simulate.hh"
#include "Interpreter.hh"
#include "Mem_Mngr.hh"
-#include "linbcg.hh"