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

Merge pull request #421 from FerhatMihoubi/master 

two commits
time-shift
Sébastien Villemot 2013-06-08 08:28:57 -07:00
parent f6f218f828 bdee6d14ec
commit 9f9c2c1ffd
4 changed files with 58 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
matlab/dr_block.m
View File

@ -301,7 +301,7 @@ for i = 1:Size;
if maximum_lead > 0 && n_fwrd > 0 if maximum_lead > 0 && n_fwrd > 0
data(i).eigval = - jacob(1 , n_pred + n - n_fwrd + 1 : n_pred + n) / jacob(1 , n_pred + n + 1 : n_pred + n + n_fwrd) ; data(i).eigval = - jacob(1 , n_pred + n - n_fwrd + 1 : n_pred + n) / jacob(1 , n_pred + n + 1 : n_pred + n + n_fwrd) ;
data(i).rank = sum(abs(data(i).eigval) > 0); data(i).rank = sum(abs(data(i).eigval) > 0);
full_rank = (abs(jacob(1,n_pred+n+1: n_pred_n+n_fwrd)) > 1e-9); full_rank = (abs(jacob(1,n_pred+n+1: n_pred+n+n_fwrd)) > 1e-9);
else else
data(i).eigval = []; data(i).eigval = [];
data(i).rank = 0; data(i).rank = 0;

10
matlab/solve_two_boundaries.m
View File

@ -68,6 +68,11 @@ g3 = [];
Blck_size=size(y_index,2); Blck_size=size(y_index,2);
correcting_factor=0.01; correcting_factor=0.01;
ilu_setup.droptol=1e-10; ilu_setup.droptol=1e-10;
ilu_setup.type = 'ilutp';
%ilu_setup.milu = 'col';
ilu_setup.milu = 'off';
ilu_setup.thresh = 1;
ilu_setup.udiag = 0;
max_resa=1e100; max_resa=1e100;
Jacobian_Size=Blck_size*(y_kmin+y_kmax_l +periods); Jacobian_Size=Blck_size*(y_kmin+y_kmax_l +periods);
g1=spalloc( Blck_size*periods, Jacobian_Size, nze*periods); g1=spalloc( Blck_size*periods, Jacobian_Size, nze*periods);
@ -244,7 +249,10 @@ while ~(cvg==1 || iter>maxit_),
elseif(stack_solve_algo==3), elseif(stack_solve_algo==3),
flag1=1; flag1=1;
while(flag1>0) while(flag1>0)
if (preconditioner == 3) if preconditioner == 2
[L1, U1]=ilu(g1a,ilu_setup);
[za,flag1] = bicgstab(g1a,b,1e-7,Blck_size*periods,L1,U1);
elseif (preconditioner == 3)
Size = Blck_size; Size = Blck_size;
gss0 = g1a(Size + 1: 2*Size,1: Size) + g1a(Size + 1: 2*Size,Size+1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size); gss0 = g1a(Size + 1: 2*Size,1: Size) + g1a(Size + 1: 2*Size,Size+1: 2*Size) + g1a(Size + 1: 2*Size,2*Size+1: 3*Size);
[L1, U1]=lu(gss0); [L1, U1]=lu(gss0);

66
mex/sources/bytecode/SparseMatrix.cc
View File

@ -183,12 +183,12 @@ dynSparseMatrix::dynSparseMatrix()
} }
dynSparseMatrix::dynSparseMatrix(const int y_size_arg, const int y_kmin_arg, const int y_kmax_arg, const bool print_it_arg, const bool steady_state_arg, const int periods_arg, dynSparseMatrix::dynSparseMatrix(const int y_size_arg, const int y_kmin_arg, const int y_kmax_arg, const bool print_it_arg, const bool steady_state_arg, const int periods_arg,
const int minimal_solving_periods_arg const int minimal_solving_periods_arg, const double slowc_arg
#ifdef CUDA #ifdef CUDA
, const int CUDA_device_arg, cublasHandle_t cublas_handle_arg, cusparseHandle_t cusparse_handle_arg, cusparseMatDescr_t descr_arg , const int CUDA_device_arg, cublasHandle_t cublas_handle_arg, cusparseHandle_t cusparse_handle_arg, cusparseMatDescr_t descr_arg
#endif #endif
): ):
Evaluate(y_size_arg, y_kmin_arg, y_kmax_arg, print_it_arg, steady_state_arg, periods_arg, minimal_solving_periods_arg) Evaluate(y_size_arg, y_kmin_arg, y_kmax_arg, print_it_arg, steady_state_arg, periods_arg, minimal_solving_periods_arg, slowc_arg)
{ {
pivotva = NULL; pivotva = NULL;
g_save_op = NULL; g_save_op = NULL;
@ -3197,6 +3197,15 @@ dynSparseMatrix::End_Solver()
End_Matlab_LU_UMFPack(); End_Matlab_LU_UMFPack();
} }
void
dynSparseMatrix::Print_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, int n)
{
int k = 0;
for (int i = 0; i < n; i++)
for (int j = Ap[i]; j < Ap[i+1]; j++)
mexPrintf("(%d, %d) %f\n", Ai[j]+1, i+1, Ax[k++]);
}
void void
dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, double *b, int n, int Size, double slowc_l, bool is_two_boundaries, int it_) dynSparseMatrix::Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, double *b, int n, int Size, double slowc_l, bool is_two_boundaries, int it_)
{ {
@ -3702,15 +3711,18 @@ dynSparseMatrix::Solve_CUDA_BiCGStab(int *Ap, int *Ai, double *Ax, int *Ap_tild,
mxSetFieldByNumber(Setup, 0, droptol, mxCreateDoubleScalar(lu_inc_tol)); mxSetFieldByNumber(Setup, 0, droptol, mxCreateDoubleScalar(lu_inc_tol));
mxSetFieldByNumber(Setup, 0, milu, mxCreateString("off")); mxSetFieldByNumber(Setup, 0, milu, mxCreateString("off"));
mxSetFieldByNumber(Setup, 0, udiag, mxCreateDoubleScalar(0)); mxSetFieldByNumber(Setup, 0, udiag, mxCreateDoubleScalar(0));
mxSetFieldByNumber(Setup, 0, thresh, mxCreateDoubleScalar(0)); mxSetFieldByNumber(Setup, 0, thresh, mxCreateDoubleScalar(1));
//mxSetFieldByNumber(Setup, 0, thresh, mxCreateDoubleScalar(1)); //mxSetFieldByNumber(Setup, 0, thresh, mxCreateDoubleScalar(1));
mxArray *lhs0[2]; mxArray *lhs0[2];
mxArray *rhs0[2]; mxArray *rhs0[2];
rhs0[0] = A_m; rhs0[0] = A_m;
rhs0[1] = Setup; rhs0[1] = Setup;
mexCallMATLAB(2, lhs0, 2, rhs0, "ilu"); ostringstream tmp;
L1 = lhs0[0]; if (mexCallMATLAB(2, lhs0, 2, rhs0, "ilu"))
U1 = lhs0[1]; {
tmp << " In BiCGStab, the incomplet LU decomposition (ilu) ahs failed.\n";
throw FatalExceptionHandling(tmp.str());
}
mxDestroyArray(Setup); mxDestroyArray(Setup);
@ -4221,7 +4233,6 @@ dynSparseMatrix::Solve_CUDA_BiCGStab(int *Ap, int *Ai, double *Ax, int *Ap_tild,
if (preconditioner == 3) if (preconditioner == 3)
{ {
double *p_tild; double *p_tild;
mexPrintf("n=%d\n",n);
cudaChk(cudaMemcpy(tmp_vect_host, p, n*sizeof(double), cudaMemcpyDeviceToHost), " in Solve_Cuda_BiCGStab, cudaMemcpy tmp_vect_host = p has failed\n"); cudaChk(cudaMemcpy(tmp_vect_host, p, n*sizeof(double), cudaMemcpyDeviceToHost), " in Solve_Cuda_BiCGStab, cudaMemcpy tmp_vect_host = p has failed\n");
/*mexPrintf("p\n"); /*mexPrintf("p\n");
@ -4716,22 +4727,24 @@ dynSparseMatrix::Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, dou
mwSize dims[1] = {(mwSize)1 }; mwSize dims[1] = {(mwSize)1 };
mxArray *Setup = mxCreateStructArray(1, dims, 5, field_names); mxArray *Setup = mxCreateStructArray(1, dims, 5, field_names);
mxSetFieldByNumber(Setup, 0, type, mxCreateString("ilutp")); mxSetFieldByNumber(Setup, 0, type, mxCreateString("ilutp"));
//mxSetFieldByNumber(Setup, 0, type, mxCreateString("nofill"));
mxSetFieldByNumber(Setup, 0, droptol, mxCreateDoubleScalar(lu_inc_tol)); mxSetFieldByNumber(Setup, 0, droptol, mxCreateDoubleScalar(lu_inc_tol));
mxSetFieldByNumber(Setup, 0, milu, mxCreateString("off")); mxSetFieldByNumber(Setup, 0, milu, mxCreateString("off"));
mxSetFieldByNumber(Setup, 0, udiag, mxCreateDoubleScalar(0)); mxSetFieldByNumber(Setup, 0, udiag, mxCreateDoubleScalar(0));
mxSetFieldByNumber(Setup, 0, thresh, mxCreateDoubleScalar(0)); mxSetFieldByNumber(Setup, 0, thresh, mxCreateDoubleScalar(1));
//mxSetFieldByNumber(Setup, 0, thresh, mxCreateDoubleScalar(1));
mxArray *lhs0[2]; mxArray *lhs0[2];
mxArray *rhs0[2]; mxArray *rhs0[2];
rhs0[0] = A_m; rhs0[0] = A_m;
rhs0[1] = Setup; rhs0[1] = Setup;
mexCallMATLAB(2, lhs0, 2, rhs0, "ilu"); if (mexCallMATLAB(2, lhs0, 2, rhs0, "ilu"))
{
ostringstream tmp;
tmp << " In BiCGStab, the incomplet LU decomposition (ilu) ahs failed.\n";
throw FatalExceptionHandling(tmp.str());
}
L1 = lhs0[0]; L1 = lhs0[0];
U1 = lhs0[1]; U1 = lhs0[1];
mxDestroyArray(Setup); mxDestroyArray(Setup);
} }
double flags = 2; double flags = 2;
mxArray *z; mxArray *z;
if (steady_state) /*Octave BicStab algorihtm involves a 0 division in case of a preconditionner equal to the LU decomposition of A matrix*/ if (steady_state) /*Octave BicStab algorihtm involves a 0 division in case of a preconditionner equal to the LU decomposition of A matrix*/
@ -6068,6 +6081,7 @@ dynSparseMatrix::Simulate_One_Boundary(int block_num, int y_size, int y_kmin, in
mxArray *b_m = NULL, *A_m = NULL, *x0_m = NULL; mxArray *b_m = NULL, *A_m = NULL, *x0_m = NULL;
SuiteSparse_long *Ap = NULL, *Ai = NULL; SuiteSparse_long *Ap = NULL, *Ai = NULL;
double *Ax = NULL, *b = NULL; double *Ax = NULL, *b = NULL;
int preconditioner = 1;
try_at_iteration = 0; try_at_iteration = 0;
@ -6138,10 +6152,12 @@ dynSparseMatrix::Simulate_One_Boundary(int block_num, int y_size, int y_kmin, in
mexPrintf("MODEL STEADY STATE: Sparse LU\n"); mexPrintf("MODEL STEADY STATE: Sparse LU\n");
break; break;
case 7: case 7:
mexPrintf("MODEL STEADY STATE: (method=GMRES)\n"); mexPrintf(preconditioner_print_out("MODEL STEADY STATE: (method=GMRES)\n", preconditioner, true).c_str());
//mexPrintf("MODEL STEADY STATE: (method=GMRES)\n");
break; break;
case 8: case 8:
mexPrintf("MODEL STEADY STATE: (method=BiCGStab)\n"); mexPrintf(preconditioner_print_out("MODEL STEADY STATE: (method=BiCGStab)\n", preconditioner, true).c_str());
//mexPrintf("MODEL STEADY STATE: (method=BiCGStab)\n");
break; break;
default: default:
mexPrintf("MODEL STEADY STATE: (method=Unknown - %d - )\n", stack_solve_algo); mexPrintf("MODEL STEADY STATE: (method=Unknown - %d - )\n", stack_solve_algo);
@ -6224,7 +6240,7 @@ dynSparseMatrix::Simulate_One_Boundary(int block_num, int y_size, int y_kmin, in
else if ((solve_algo == 7 && steady_state) || (stack_solve_algo == 2 && !steady_state)) else if ((solve_algo == 7 && steady_state) || (stack_solve_algo == 2 && !steady_state))
Solve_Matlab_GMRES(A_m, b_m, size, slowc, block_num, false, it_, x0_m); Solve_Matlab_GMRES(A_m, b_m, size, slowc, block_num, false, it_, x0_m);
else if ((solve_algo == 8 && steady_state) || (stack_solve_algo == 3 && !steady_state)) else if ((solve_algo == 8 && steady_state) || (stack_solve_algo == 3 && !steady_state))
Solve_Matlab_BiCGStab(A_m, b_m, size, slowc, block_num, false, it_, x0_m, 1); Solve_Matlab_BiCGStab(A_m, b_m, size, slowc, block_num, false, it_, x0_m, preconditioner);
else if ((solve_algo == 6 && steady_state) || ((stack_solve_algo == 0 || stack_solve_algo == 1 || stack_solve_algo == 4) && !steady_state)) else if ((solve_algo == 6 && steady_state) || ((stack_solve_algo == 0 || stack_solve_algo == 1 || stack_solve_algo == 4) && !steady_state))
Solve_LU_UMFPack(Ap, Ai, Ax, b, size, size, slowc, true, 0); Solve_LU_UMFPack(Ap, Ai, Ax, b, size, size, slowc, true, 0);
} }
@ -6339,20 +6355,26 @@ dynSparseMatrix::Simulate_Newton_One_Boundary(const bool forward)
} }
string string
dynSparseMatrix::preconditioner_print_out(string s, int preconditioner) dynSparseMatrix::preconditioner_print_out(string s, int preconditioner, bool ss)
{ {
int n = s.length(); int n = s.length();
string tmp = ", preconditioner="; string tmp = ", preconditioner=";
switch(preconditioner) switch(preconditioner)
{ {
case 0: case 0:
tmp.append("Jacobi on dynamic jacobian"); if (ss)
tmp.append("Jacobi on static jacobian");
else
tmp.append("Jacobi on dynamic jacobian");
break; break;
case 1: case 1:
tmp.append("incomplet lu0 on dynamic jacobian"); if (ss)
tmp.append("incomplet lutp on static jacobian");
else
tmp.append("incomplet lu0 on dynamic jacobian");
break; break;
case 2: case 2:
tmp.append("incomplet lut on dynamic jacobian"); tmp.append("incomplet lutp on dynamic jacobian");
break; break;
case 3: case 3:
tmp.append("lu on static jacobian"); tmp.append("lu on static jacobian");
@ -6526,10 +6548,10 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin
mexPrintf("MODEL SIMULATION: (method=Relaxation)\n"); mexPrintf("MODEL SIMULATION: (method=Relaxation)\n");
break; break;
case 2: case 2:
mexPrintf(preconditioner_print_out("MODEL SIMULATION: (method=GMRES)\n", preconditioner).c_str()); mexPrintf(preconditioner_print_out("MODEL SIMULATION: (method=GMRES)\n", preconditioner, false).c_str());
break; break;
case 3: case 3:
mexPrintf(preconditioner_print_out("MODEL SIMULATION: (method=BiCGStab)\n", preconditioner).c_str()); mexPrintf(preconditioner_print_out("MODEL SIMULATION: (method=BiCGStab)\n", preconditioner, false).c_str());
break; break;
case 4: case 4:
mexPrintf("MODEL SIMULATION: (method=Sparse LU & optimal path length)\n"); mexPrintf("MODEL SIMULATION: (method=Sparse LU & optimal path length)\n");
@ -6538,7 +6560,7 @@ dynSparseMatrix::Simulate_Newton_Two_Boundaries(int blck, int y_size, int y_kmin
mexPrintf("MODEL SIMULATION: (method=ByteCode own solver)\n"); mexPrintf("MODEL SIMULATION: (method=ByteCode own solver)\n");
break; break;
case 7: case 7:
mexPrintf(preconditioner_print_out("MODEL SIMULATION: (method=GPU BiCGStab)\n", preconditioner).c_str()); mexPrintf(preconditioner_print_out("MODEL SIMULATION: (method=GPU BiCGStab)\n", preconditioner, false).c_str());
break; break;
default: default:
mexPrintf("MODEL SIMULATION: (method=Unknown - %d - )\n", stack_solve_algo); mexPrintf("MODEL SIMULATION: (method=Unknown - %d - )\n", stack_solve_algo);

7
mex/sources/bytecode/SparseMatrix.hh
View File

@ -108,7 +108,7 @@ public:
typedef int64_t SuiteSparse_long; typedef int64_t SuiteSparse_long;
#endif #endif
dynSparseMatrix(); dynSparseMatrix();
dynSparseMatrix(const int y_size_arg, const int y_kmin_arg, const int y_kmax_arg, const bool print_it_arg, const bool steady_state_arg, const int periods_arg, const int minimal_solving_periods_arg dynSparseMatrix(const int y_size_arg, const int y_kmin_arg, const int y_kmax_arg, const bool print_it_arg, const bool steady_state_arg, const int periods_arg, const int minimal_solving_periods_arg, const double slowc_arg
#ifdef CUDA #ifdef CUDA
,const int CUDA_device_arg, cublasHandle_t cublas_handle_arg, cusparseHandle_t cusparse_handle_arg, cusparseMatDescr_t descr_arg ,const int CUDA_device_arg, cublasHandle_t cublas_handle_arg, cusparseHandle_t cusparse_handle_arg, cusparseMatDescr_t descr_arg
#endif #endif
@ -141,6 +141,7 @@ private:
bool Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, int it_); bool Solve_ByteCode_Sparse_GaussianElimination(int Size, int blck, int it_);
void Solve_Matlab_Relaxation(mxArray *A_m, mxArray *b_m, unsigned int Size, double slowc_l, bool is_two_boundaries, int it_); void Solve_Matlab_Relaxation(mxArray *A_m, mxArray *b_m, unsigned int Size, double slowc_l, bool is_two_boundaries, int it_);
void Solve_Matlab_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double slowc_l, bool is_two_boundaries, int it_); void Solve_Matlab_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double slowc_l, bool is_two_boundaries, int it_);
void Print_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, int n);
void Solve_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double slowc_l, bool is_two_boundaries, int it_); void Solve_LU_UMFPack(mxArray *A_m, mxArray *b_m, int Size, double slowc_l, bool is_two_boundaries, int it_);
void Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, double *b, int n, int Size, double slowc_l, bool is_two_boundaries, int it_); void Solve_LU_UMFPack(SuiteSparse_long *Ap, SuiteSparse_long *Ai, double *Ax, double *b, int n, int Size, double slowc_l, bool is_two_boundaries, int it_);
void End_Matlab_LU_UMFPack(); void End_Matlab_LU_UMFPack();
@ -157,7 +158,7 @@ private:
bool Simulate_One_Boundary(int blck, int y_size, int y_kmin, int y_kmax, int Size, bool cvg); bool Simulate_One_Boundary(int blck, int y_size, int y_kmin, int y_kmax, int Size, bool cvg);
bool solve_linear(const int block_num, const int y_size, const int y_kmin, const int y_kmax, const int size, const int iter); bool solve_linear(const int block_num, const int y_size, const int y_kmin, const int y_kmax, const int size, const int iter);
void solve_non_linear(const int block_num, const int y_size, const int y_kmin, const int y_kmax, const int size); void solve_non_linear(const int block_num, const int y_size, const int y_kmin, const int y_kmax, const int size);
string preconditioner_print_out(string s, int preconditioner); string preconditioner_print_out(string s, int preconditioner, bool ss);
bool compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4, int Size bool compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4, int Size
#ifdef PROFILER #ifdef PROFILER
, long int *ndiv, long int *nsub , long int *ndiv, long int *nsub
@ -237,7 +238,7 @@ protected:
int u_count_alloc, u_count_alloc_save; int u_count_alloc, u_count_alloc_save;
vector<double *> jac; vector<double *> jac;
double *jcb; double *jcb;
double slowc, slowc_save, prev_slowc_save, markowitz_c; double slowc_save, prev_slowc_save, markowitz_c;
int y_decal; int y_decal;
int *index_equa; int *index_equa;
int u_count, tbreak_g; int u_count, tbreak_g;