SparseMatrix.hh 6.22 KB
Newer Older
sebastien's avatar
sebastien committed
1
/*
sebastien's avatar
trunk:    
sebastien committed
2
 * Copyright (C) 2007-2009 Dynare Team
sebastien's avatar
sebastien committed
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
 *
 * 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/>.
 */

ferhat's avatar
ferhat committed
20
21
22
23
24
#ifndef SPARSEMATRIX_HH_INCLUDED
#define SPARSEMATRIX_HH_INCLUDED

#include <fstream>
#include <stack>
25
#include <cmath>
ferhat's avatar
ferhat committed
26
#include <map>
27
#include <ctime>
ferhat's avatar
ferhat committed
28
29
30
31
32
#include "Mem_Mngr.hh"
#define NEW_ALLOC
#define MARKOVITZ

using namespace std;
33
34
35

#ifdef _MSC_VER
# include <limits>
36
37
38
39
#include <boost/math/special_functions/erf.hpp>
#define M_SQRT2 1.4142135623730950488016887242097
#define M_PI 3.1415926535897932384626433832795
#define erf(x) boost::math::erf(x)
40
41
42
extern unsigned long _nan[2];
extern double NAN;

sebastien's avatar
sebastien committed
43
44
inline bool
isnan(double value)
45
46
47
48
{
  return _isnan(value);
}

sebastien's avatar
sebastien committed
49
50
inline bool
isinf(double value)
51
{
sebastien's avatar
sebastien committed
52
53
  return (std::numeric_limits<double>::has_infinity
          && value == std::numeric_limits<double>::infinity());
54
55
56
}

template<typename T>
sebastien's avatar
sebastien committed
57
58
inline T
asinh(T x)
59
60
61
62
63
{
  return log(x+sqrt(x*x+1));
}

template<typename T>
sebastien's avatar
sebastien committed
64
65
inline T
acosh(T x)
66
{
sebastien's avatar
sebastien committed
67
  if (!(x >= 1.0))
68
69
70
71
72
    return sqrt(-1.0);
  return log(x+sqrt(x*x-1.0));
}

template<typename T>
sebastien's avatar
sebastien committed
73
74
inline T
atanh(T x)
75
{
sebastien's avatar
sebastien committed
76
  if (!(x > -1.0 && x < 1.0))
77
78
79
80
81
    return sqrt(-1.0);
  return log((1.0+x)/(1.0-x))/2.0;
}

#endif
ferhat's avatar
ferhat committed
82

83
struct t_save_op_s
ferhat's avatar
ferhat committed
84
85
{
  short int lag, operat;
86
  int first, second;
ferhat's avatar
ferhat committed
87
88
};

sebastien's avatar
sebastien committed
89
90
91
92
93
94
95
96
97
98
99
const int IFLD  = 0;
const int IFDIV = 1;
const int IFLESS = 2;
const int IFSUB = 3;
const int IFLDZ = 4;
const int IFMUL = 5;
const int IFSTP = 6;
const int IFADD = 7;
const double eps = 1e-10;
const double very_big = 1e24;
const int alt_symbolic_count_max = 1;
100
const double mem_increasing_factor = 1.1;
ferhat's avatar
ferhat committed
101
102

class SparseMatrix
sebastien's avatar
sebastien committed
103
104
105
{
public:
  SparseMatrix();
106
107
  int simulate_NG1(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, bool cvg, int &iter, int minimal_solving_periods, int Block_number, int stack_solve_algo);
  bool simulate_NG(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, bool print_it, bool cvg, int &iter, bool steady_state, int Block_number, int solve_algo);
sebastien's avatar
sebastien committed
108
109
  void Direct_Simulate(int blck, int y_size, int it_, int y_kmin, int y_kmax, int Size, int periods, bool print_it, int iter);
  void fixe_u(double **u, int u_count_int, int max_lag_plus_max_lead_plus_1);
110
  void Read_SparseMatrix(string file_name, const int Size, int periods, int y_kmin, int y_kmax, bool steady_state, bool two_boundaries, int stack_solve_algo, int solve_algo);
sebastien's avatar
sebastien committed
111
  void Read_file(string file_name, int periods, int u_size1, int y_size, int y_kmin, int y_kmax, int &nb_endo, int &u_count, int &u_count_init, double *u);
112
  double g0, gp0, glambda2, try_at_iteration;
113

sebastien's avatar
sebastien committed
114
private:
115
116
117
  void Init_GE(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM);
  void Init_Matlab_Sparse(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m);
  void Init_Matlab_Sparse_Simple(int Size, map<pair<pair<int, int>, int>, int> &IM, mxArray *A_m, mxArray *b_m);
sebastien's avatar
sebastien committed
118
  void Simple_Init(int it_, int y_kmin, int y_kmax, int Size, std::map<std::pair<std::pair<int, int>, int>, int> &IM);
119
120
121
122
  void End_GE(int Size);
  void Solve_Matlab_LU_UMFPack(mxArray* A_m, mxArray* b_m, int Size, double slowc_l);
  void Solve_Matlab_GMRES(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block);
  void Solve_Matlab_BiCGStab(mxArray* A_m, mxArray* b_m, int Size, double slowc, int block);
sebastien's avatar
sebastien committed
123
  bool compare(int *save_op, int *save_opa, int *save_opaa, int beg_t, int periods, long int nop4,  int Size
ferhat's avatar
ferhat committed
124
#ifdef PROFILER
sebastien's avatar
sebastien committed
125
               , long int *ndiv, long int *nsub
ferhat's avatar
ferhat committed
126
#endif
sebastien's avatar
sebastien committed
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
               );
  void Insert(const int r, const int c, const int u_index, const int lag_index);
  void Delete(const int r, const int c);
  int At_Row(int r, NonZeroElem **first);
  int At_Pos(int r, int c, NonZeroElem **first);
  int At_Col(int c, NonZeroElem **first);
  int At_Col(int c, int lag, NonZeroElem **first);
  int NRow(int r);
  int NCol(int c);
  int Union_Row(int row1, int row2);
  void Print(int Size, int *b);
  int Get_u();
  void Delete_u(int pos);
  void Clear_u();
  void Print_u();
  void CheckIt(int y_size, int y_kmin, int y_kmax, int Size, int periods, int iter);
  void Check_the_Solution(int periods, int y_kmin, int y_kmax, int Size, double *u, int *pivot, int *b);
  int complete(int beg_t, int Size, int periods, int *b);
  double bksub(int tbreak, int last_period, int Size, double slowc_l
ferhat's avatar
ferhat committed
146
#ifdef PROFILER
sebastien's avatar
sebastien committed
147
               , long int *nmul
ferhat's avatar
ferhat committed
148
#endif
sebastien's avatar
sebastien committed
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
               );
  double simple_bksub(int it_, int Size, double slowc_l);
  stack<double> Stack;
  int nb_prologue_table_u, nb_first_table_u, nb_middle_table_u, nb_last_table_u;
  int nb_prologue_table_y, nb_first_table_y, nb_middle_table_y, nb_last_table_y;
  int middle_count_loop;
  char type;
  fstream SaveCode;
  string filename;
  int max_u, min_u;
  clock_t time00;

  Mem_Mngr mem_mngr;
  vector<int> u_liste;
  map<pair<int, int>, NonZeroElem *> Mapped_Array;
  int *NbNZRow, *NbNZCol;
  NonZeroElem **FNZE_R, **FNZE_C;
  int nb_endo, u_count_init;

  int *pivot, *pivotk, *pivot_save;
  double *pivotv, *pivotva;
  int *b;
  bool *line_done;
  bool symbolic, alt_symbolic;
  int alt_symbolic_count;
  int *g_save_op;
  int first_count_loop;
  int g_nop_all;
  double markowitz_c_s;
  double res1a;
  long int nop_all, nop1, nop2;
  map<pair<pair<int, int>, int>, int> IM_i;
ferhat's avatar
ferhat committed
181
protected:
182
  int u_count_alloc, u_count_alloc_save;
sebastien's avatar
sebastien committed
183
  double *u, *y, *ya;
184
185
  vector<double*> jac;
  double *jcb;
sebastien's avatar
sebastien committed
186
  double res1, res2, max_res, max_res_idx;
187
  double slowc, slowc_save, prev_slowc_save, markowitz_c;
sebastien's avatar
sebastien committed
188
189
190
191
192
193
194
195
  int y_kmin, y_kmax, y_size, periods, y_decal;
  int  *index_vara, *index_equa;
  int u_count, tbreak_g;
  int iter;
  double *direction;
  int start_compare;
  int restart;
  bool error_not_printed;
196
  double g_lambda1, g_lambda2, gp_0;
197
198
  double lu_inc_tol;
  bool reduced;
sebastien's avatar
sebastien committed
199
};
ferhat's avatar
ferhat committed
200
201

#endif