SparseMatrix.cc 64.7 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/>.
 */

20
#include <cstring>
21
#include <ctime>
22
#include <sstream>
ferhat's avatar
ferhat committed
23
24
#include "SparseMatrix.hh"

25
26
27
28
#ifdef _MSC_VER
unsigned long _nan[2] = { 0xffffffff, 0x7fffffff };
double NAN = *((double *) _nan);
#endif
29

ferhat's avatar
ferhat committed
30
31
SparseMatrix::SparseMatrix()
{
32
33
34
  pivotva = NULL;
  g_save_op = NULL;
  g_nop_all = 0;
ferhat's avatar
ferhat committed
35
  mem_mngr.init_Mem();
36
37
38
39
40
41
42
43
  symbolic = true;
  alt_symbolic = false;
  alt_symbolic_count = 0;
  max_u = 0;
  min_u = 0x7FFFFFFF;
  res1a = 9.0e60;
  tbreak_g = 0;
  start_compare = 0;
44
  restart = 0;
45
  IM_i.clear();
ferhat's avatar
ferhat committed
46
47
}

48
49
int
SparseMatrix::NRow(int r)
ferhat's avatar
ferhat committed
50
51
52
53
{
  return NbNZRow[r];
}

54
55
int
SparseMatrix::NCol(int c)
ferhat's avatar
ferhat committed
56
57
58
59
{
  return NbNZCol[c];
}

60
61
int
SparseMatrix::At_Row(int r, NonZeroElem **first)
ferhat's avatar
ferhat committed
62
{
63
  (*first) = FNZE_R[r];
ferhat's avatar
ferhat committed
64
65
66
67
68
69
70
  return NbNZRow[r];
}

int
SparseMatrix::Union_Row(int row1, int row2)
{
  NonZeroElem *first1, *first2;
71
72
73
74
  int n1 = At_Row(row1, &first1);
  int n2 = At_Row(row2, &first2);
  int i1 = 0, i2 = 0, nb_elem = 0;
  while (i1 < n1 && i2 < n2)
ferhat's avatar
ferhat committed
75
    {
76
      if (first1->c_index == first2->c_index)
ferhat's avatar
ferhat committed
77
78
79
80
        {
          nb_elem++;
          i1++;
          i2++;
81
82
          first1 = first1->NZE_R_N;
          first2 = first2->NZE_R_N;
ferhat's avatar
ferhat committed
83
        }
84
      else if (first1->c_index < first2->c_index)
ferhat's avatar
ferhat committed
85
86
87
        {
          nb_elem++;
          i1++;
88
          first1 = first1->NZE_R_N;
ferhat's avatar
ferhat committed
89
90
91
92
93
        }
      else
        {
          nb_elem++;
          i2++;
94
          first2 = first2->NZE_R_N;
ferhat's avatar
ferhat committed
95
96
97
98
99
100
101
102
        }
    }
  return nb_elem;
}

int
SparseMatrix::At_Pos(int r, int c, NonZeroElem **first)
{
103
104
105
  (*first) = FNZE_R[r];
  while ((*first)->c_index != c)
    (*first) = (*first)->NZE_R_N;
ferhat's avatar
ferhat committed
106
107
108
  return NbNZRow[r];
}

109
110
int
SparseMatrix::At_Col(int c, NonZeroElem **first)
ferhat's avatar
ferhat committed
111
{
112
  (*first) = FNZE_C[c];
ferhat's avatar
ferhat committed
113
114
115
  return NbNZCol[c];
}

116
117
int
SparseMatrix::At_Col(int c, int lag, NonZeroElem **first)
ferhat's avatar
ferhat committed
118
{
119
120
121
122
  (*first) = FNZE_C[c];
  int i = 0;
  while ((*first)->lag_index != lag && (*first))
    (*first) = (*first)->NZE_C_N;
ferhat's avatar
ferhat committed
123
124
  if ((*first))
    {
125
      NonZeroElem *firsta = (*first);
ferhat's avatar
ferhat committed
126
127
128
129
      if (!firsta->NZE_C_N)
        i++;
      else
        {
130
          while (firsta->lag_index == lag && firsta->NZE_C_N)
ferhat's avatar
ferhat committed
131
            {
132
              firsta = firsta->NZE_C_N;
ferhat's avatar
ferhat committed
133
134
              i++;
            }
135
          if (firsta->lag_index == lag) i++;
ferhat's avatar
ferhat committed
136
137
138
139
140
        }
    }
  return i;
}

141
142
void
SparseMatrix::Delete(const int r, const int c)
ferhat's avatar
ferhat committed
143
{
144
145
146
147
148
149
150
151
152
153
154
	NonZeroElem *first = FNZE_R[r], *firsta = NULL;

  while (first->c_index != c)
    {
      firsta = first;
      first = first->NZE_R_N;
    }
  if (firsta != NULL)
    firsta->NZE_R_N = first->NZE_R_N;
  if (first == FNZE_R[r])
    FNZE_R[r] = first->NZE_R_N;
ferhat's avatar
ferhat committed
155
  NbNZRow[r]--;
156
157
158
159

  first = FNZE_C[c];
  firsta = NULL;
  while (first->r_index != r)
ferhat's avatar
ferhat committed
160
    {
161
162
      firsta = first;
      first = first->NZE_C_N;
ferhat's avatar
ferhat committed
163
    }
164
165
166
167
168
169

  if (firsta != NULL)
    firsta->NZE_C_N = first->NZE_C_N;
  if (first == FNZE_C[c])
		FNZE_C[c] = first->NZE_C_N;

ferhat's avatar
ferhat committed
170
171
172
173
174
  u_liste.push_back(first->u_index);
  mem_mngr.mxFree_NZE(first);
  NbNZCol[c]--;
}

175
176
void
SparseMatrix::Print(int Size, int *b)
ferhat's avatar
ferhat committed
177
{
178
  int a, i, j, k, l;
ferhat's avatar
ferhat committed
179
  mexPrintf("   ");
180
181
  for (k = 0; k < Size*periods; k++)
    mexPrintf("%-2d ", k);
ferhat's avatar
ferhat committed
182
  mexPrintf("    |    ");
183
184
  for (k = 0; k < Size*periods; k++)
    mexPrintf("%8d", k);
ferhat's avatar
ferhat committed
185
  mexPrintf("\n");
186
  for (i = 0; i < Size*periods; i++)
ferhat's avatar
ferhat committed
187
    {
188
189
190
191
192
      NonZeroElem *first = FNZE_R[i];
      j = NbNZRow[i];
      mexPrintf("%-2d ", i);
      a = 0;
      for (k = 0; k < j; k++)
ferhat's avatar
ferhat committed
193
        {
194
          for (l = 0; l < (first->c_index-a); l++)
ferhat's avatar
ferhat committed
195
            mexPrintf("   ");
196
197
198
          mexPrintf("%-2d ", first->u_index);
          a = first->c_index+1;
          first = first->NZE_R_N;
ferhat's avatar
ferhat committed
199
        }
200
      for (k = a; k < Size*periods; k++)
ferhat's avatar
ferhat committed
201
        mexPrintf("   ");
202
      mexPrintf("%-2d ", b[i]);
ferhat's avatar
ferhat committed
203

204
205
206
207
208
      first = FNZE_R[i];
      j = NbNZRow[i];
      mexPrintf(" | %-2d ", i);
      a = 0;
      for (k = 0; k < j; k++)
ferhat's avatar
ferhat committed
209
        {
210
          for (l = 0; l < (first->c_index-a); l++)
ferhat's avatar
ferhat committed
211
            mexPrintf("        ");
212
213
214
          mexPrintf("%8.4f", double (u[first->u_index]));
          a = first->c_index+1;
          first = first->NZE_R_N;
ferhat's avatar
ferhat committed
215
        }
216
      for (k = a; k < Size*periods; k++)
ferhat's avatar
ferhat committed
217
        mexPrintf("        ");
218
      mexPrintf("%8.4f", double (u[b[i]]));
ferhat's avatar
ferhat committed
219
220
221
222
      mexPrintf("\n");
    }
}

223
224
void
SparseMatrix::Insert(const int r, const int c, const int u_index, const int lag_index)
ferhat's avatar
ferhat committed
225
{
226
227
228
229
230
  NonZeroElem *firstn, *first, *firsta, *a;
  firstn = mem_mngr.mxMalloc_NZE();
  first = FNZE_R[r];
  firsta = NULL;
  while (first->c_index < c && (a = first->NZE_R_N))
231
    {
232
233
      firsta = first;
      first = a;
ferhat's avatar
ferhat committed
234
    }
235
236
237
238
239
  firstn->u_index = u_index;
  firstn->r_index = r;
  firstn->c_index = c;
  firstn->lag_index = lag_index;
  if (first->c_index > c)
240
    {
241
242
243
244
245
      if (first == FNZE_R[r])
        FNZE_R[r] = firstn;
      if (firsta != NULL)
        firsta->NZE_R_N = firstn;
      firstn->NZE_R_N = first;
ferhat's avatar
ferhat committed
246
    }
247
  else
248
    {
249
250
      first->NZE_R_N = firstn;
      firstn->NZE_R_N = NULL;
251
252
    }
  NbNZRow[r]++;
253
254
255
  first = FNZE_C[c];
  firsta = NULL;
  while (first->r_index < r && (a = first->NZE_C_N))
ferhat's avatar
ferhat committed
256
    {
257
258
      firsta = first;
      first = a;
ferhat's avatar
ferhat committed
259
    }
260
  if (first->r_index > r)
ferhat's avatar
ferhat committed
261
    {
262
263
264
265
266
      if (first == FNZE_C[c])
        FNZE_C[c] = firstn;
      if (firsta != NULL)
        firsta->NZE_C_N = firstn;
      firstn->NZE_C_N = first;
ferhat's avatar
ferhat committed
267
    }
268
  else
ferhat's avatar
ferhat committed
269
    {
270
271
      first->NZE_C_N = firstn;
      firstn->NZE_C_N = NULL;
ferhat's avatar
ferhat committed
272
    }
273

ferhat's avatar
ferhat committed
274
275
276
  NbNZCol[c]++;
}

277
void
278
SparseMatrix::Read_SparseMatrix(string file_name, const int Size, int periods, int y_kmin, int y_kmax, bool steady_state, bool two_boundaries)
ferhat's avatar
ferhat committed
279
{
280
281
  unsigned int eq, var;
  int i, j, lag;
282
  filename = file_name;
ferhat's avatar
ferhat committed
283
284
285
  mem_mngr.fixe_file_name(file_name);
  if (!SaveCode.is_open())
    {
286
287
288
289
    	if(steady_state)
        SaveCode.open((file_name + "_static.bin").c_str(), ios::in | ios::binary);
			else
			  SaveCode.open((file_name + "_dynamic.bin").c_str(), ios::in | ios::binary);
ferhat's avatar
ferhat committed
290
291
      if (!SaveCode.is_open())
        {
292
293
294
295
        	if(steady_state)
            mexPrintf("Error : Can't open file \"%s\" for reading\n", (file_name + "_static.bin").c_str());
					else
					  mexPrintf("Error : Can't open file \"%s\" for reading\n", (file_name + "_dynamic.bin").c_str());
ferhat's avatar
ferhat committed
296
297
298
299
300
          mexEvalString("st=fclose('all');clear all;");
          mexErrMsgTxt("Exit from Dynare");
        }
    }
  IM_i.clear();
301
302
303
304
305
306
307
308
309
310
311
	if(two_boundaries)
	  {
	  	for (i = 0; i < u_count_init-Size; i++)
        {
          SaveCode.read(reinterpret_cast<char *>(&eq), sizeof(eq));
          SaveCode.read(reinterpret_cast<char *>(&var), sizeof(var));
          SaveCode.read(reinterpret_cast<char *>(&lag), sizeof(lag));
          SaveCode.read(reinterpret_cast<char *>(&j), sizeof(j));
          IM_i[make_pair(make_pair(eq, var), lag)] = j;
        }
      for (j=0;j<Size;j++)
312
        IM_i[make_pair(make_pair(j, Size*(periods+y_kmax)), 0)] = j;
313
314
315
316
317
318
319
320
321
322
323
324
	  }
	else
	  {
	  	for (i = 0; i < u_count_init; i++)
        {
          SaveCode.read(reinterpret_cast<char *>(&eq), sizeof(eq));
          SaveCode.read(reinterpret_cast<char *>(&var), sizeof(var));
          SaveCode.read(reinterpret_cast<char *>(&lag), sizeof(lag));
          SaveCode.read(reinterpret_cast<char *>(&j), sizeof(j));
          IM_i[make_pair(make_pair(eq, var), lag)] = j;
        }
	  }
325
326
327
328
  index_vara = (int *) mxMalloc(Size*(periods+y_kmin+y_kmax)*sizeof(int));
  for (j = 0; j < Size; j++)
    SaveCode.read(reinterpret_cast<char *>(&index_vara[j]), sizeof(*index_vara));
  if (periods+y_kmin+y_kmax > 1)
ferhat's avatar
ferhat committed
329
    {
330
      for (i = 1; i < periods+y_kmin+y_kmax; i++)
ferhat's avatar
ferhat committed
331
        {
332
333
334
335
336
          for (j = 0; j < Size; j++)
            {
              index_vara[j+Size*i] = index_vara[j+Size*(i-1)]+y_size;
            }
        }
ferhat's avatar
ferhat committed
337
    }
338
339
  index_equa = (int *) mxMalloc(Size*sizeof(int));
  for (j = 0; j < Size; j++)
ferhat's avatar
ferhat committed
340
341
342
343
344
    {
      SaveCode.read(reinterpret_cast<char *>(&index_equa[j]), sizeof(*index_equa));
    }
}

345
346
void
SparseMatrix::Simple_Init(int it_, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM)
347
348
{
  int i, eq, var, lag;
349
350
351
352
353
354
355
356
357
  map<pair<pair<int, int>, int>, int>::iterator it4;
  NonZeroElem *first;
  pivot = (int *) mxMalloc(Size*sizeof(int));
  pivot_save = (int *) mxMalloc(Size*sizeof(int));
  pivotk = (int *) mxMalloc(Size*sizeof(int));
  pivotv = (double *) mxMalloc(Size*sizeof(double));
  pivotva = (double *) mxMalloc(Size*sizeof(double));
  b = (int *) mxMalloc(Size*sizeof(int));
  line_done = (bool *) mxMalloc(Size*sizeof(bool));
358

359
  mem_mngr.init_CHUNK_BLCK_SIZE(u_count);
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
  g_save_op = NULL;
  g_nop_all = 0;
  i = Size*sizeof(NonZeroElem *);
  FNZE_R = (NonZeroElem **) mxMalloc(i);
  FNZE_C = (NonZeroElem **) mxMalloc(i);
  NonZeroElem **temp_NZE_R = (NonZeroElem **) mxMalloc(i);
  NonZeroElem **temp_NZE_C = (NonZeroElem **) mxMalloc(i);
  i = Size*sizeof(int);
  NbNZRow = (int *) mxMalloc(i);
  NbNZCol = (int *) mxMalloc(i);
  i = Size*sizeof(*b);
  it4 = IM.begin();
  eq = -1;
  //#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
  for (i = 0; i < Size; i++)
375
    {
376
377
378
379
380
381
382
383
      b[i] = 0;
      line_done[i] = 0;
      FNZE_C[i] = 0;
      FNZE_R[i] = 0;
      temp_NZE_C[i] = 0;
      temp_NZE_R[i] = 0;
      NbNZRow[i] = 0;
      NbNZCol[i] = 0;
384
    }
385
386
  int u_count1 = Size;
  while (it4 != IM.end())
387
    {
388
389
390
391
      var = it4->first.first.second;
      eq = it4->first.first.first;
      lag = it4->first.second;
      if (lag == 0)   /*Build the index for sparse matrix containing the jacobian : u*/
392
393
394
        {
          NbNZRow[eq]++;
          NbNZCol[var]++;
395
396
397
          first = mem_mngr.mxMalloc_NZE();
          first->NZE_C_N = NULL;
          first->NZE_R_N = NULL;
398
          first->u_index = u_count1;
399
400
401
402
403
404
405
406
407
408
409
410
411
          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;
412
413
414
415
          u_count1++;
        }
      it4++;
    }
416
417
418
  //#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
  for (i = 0; i < Size; i++)
		b[i] = i;
419
420
  mxFree(temp_NZE_R);
  mxFree(temp_NZE_C);
421
  u_count = u_count1;
422
423
}

424
425
void
SparseMatrix::Init(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM)
ferhat's avatar
ferhat committed
426
{
427
428
429
430
431
432
433
434
435
436
437
  int t, i, eq, var, lag, ti_y_kmin, ti_y_kmax;
  double tmp_b = 0.0;
  map<pair<pair<int, int>, int>, int>::iterator it4;
  NonZeroElem *first;
  pivot = (int *) mxMalloc(Size*periods*sizeof(int));
  pivot_save = (int *) mxMalloc(Size*periods*sizeof(int));
  pivotk = (int *) mxMalloc(Size*periods*sizeof(int));
  pivotv = (double *) mxMalloc(Size*periods*sizeof(double));
  pivotva = (double *) mxMalloc(Size*periods*sizeof(double));
  b = (int *) mxMalloc(Size*periods*sizeof(int));
  line_done = (bool *) mxMalloc(Size*periods*sizeof(bool));
ferhat's avatar
ferhat committed
438
  mem_mngr.init_CHUNK_BLCK_SIZE(u_count);
439
440
441
442
443
444
445
446
447
448
449
450
451
  g_save_op = NULL;
  g_nop_all = 0;
  i = (periods+y_kmax+1)*Size*sizeof(NonZeroElem*);
  FNZE_R = (NonZeroElem **) mxMalloc(i);
  FNZE_C = (NonZeroElem **) mxMalloc(i);
  NonZeroElem **temp_NZE_R = (NonZeroElem **) mxMalloc(i);
  NonZeroElem **temp_NZE_C = (NonZeroElem **) mxMalloc(i);
  i = (periods+y_kmax+1)*Size*sizeof(int);
  NbNZRow = (int *) mxMalloc(i);
  NbNZCol = (int *) mxMalloc(i);

  //#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
  for (i = 0; i < periods*Size; i++)
452
    {
453
454
      b[i] = 0;
      line_done[i] = 0;
455
    }
456
457
  //#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
  for (i = 0; i < (periods+y_kmax+1)*Size; i++)
458
    {
459
460
461
462
463
464
      FNZE_C[i] = 0;
      FNZE_R[i] = 0;
      temp_NZE_C[i] = NULL;
      temp_NZE_R[i] = NULL;
      NbNZRow[i] = 0;
      NbNZCol[i] = 0;
465
466
    }

467
468
  //#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) ordered private(it4, ti_y_kmin, ti_y_kmax, eq, var, lag) schedule(dynamic)
  for (t = 0; t < periods; t++)
ferhat's avatar
ferhat committed
469
    {
470
471
472
473
474
475
      ti_y_kmin = -min(t, y_kmin);
      ti_y_kmax = min(periods-(t+1), y_kmax);
      it4 = IM.begin();
      eq = -1;
      //#pragma omp ordered
      while (it4 != IM.end())
ferhat's avatar
ferhat committed
476
        {
477
478
479
480
481
482
          var = it4->first.first.second;
          if (eq != it4->first.first.first+Size*t)
            tmp_b = 0;
          eq = it4->first.first.first+Size*t;
          lag = it4->first.second;
          if (var < (periods+y_kmax)*Size)
ferhat's avatar
ferhat committed
483
            {
484
485
              lag = it4->first.second;
              if (lag <= ti_y_kmax && lag >= ti_y_kmin)   /*Build the index for sparse matrix containing the jacobian : u*/
ferhat's avatar
ferhat committed
486
                {
487
                  var += Size*t;
ferhat's avatar
ferhat committed
488
489
                  NbNZRow[eq]++;
                  NbNZCol[var]++;
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
                  first = mem_mngr.mxMalloc_NZE();
                  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;
ferhat's avatar
ferhat committed
507
                }
508
              else       /*Build the additive terms ooutside the simulation periods related to the first lags and the last leads...*/
ferhat's avatar
ferhat committed
509
                {
510
511
512
513
514
515
516
                  if (lag < ti_y_kmin)
                    {
                      tmp_b += u[it4->second+u_count_init*t]*y[index_vara[var+Size*(y_kmin+t)]];
                    }
                  else
                    {
                      tmp_b += u[it4->second+u_count_init*t]*y[index_vara[var+Size*(y_kmin+t)]];
517

518
                    }
ferhat's avatar
ferhat committed
519
520
                }
            }
521
          else           /* ...and store it in the u vector*/
ferhat's avatar
ferhat committed
522
            {
523
524
              b[eq] = it4->second+u_count_init*t;
              u[b[eq]] += tmp_b;
525
              tmp_b = 0;
ferhat's avatar
ferhat committed
526
527
528
529
530
531
532
533
            }
          it4++;
        }
    }
  mxFree(temp_NZE_R);
  mxFree(temp_NZE_C);
}

534
535
void
SparseMatrix::ShortInit(int periods, int y_kmin, int y_kmax, int Size, map<pair<pair<int, int>, int>, int> &IM)
ferhat's avatar
ferhat committed
536
{
537
  int t, eq, var, lag, ti_y_kmin, ti_y_kmax;
538
539
  double tmp_b = 0.0;
  map<pair<pair<int, int>, int>, int>::iterator it4;
540
  //#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) ordered private(it4, ti_y_kmin, ti_y_kmax, eq, var, lag, tmp_b) schedule(dynamic)
541
  for (t = 0; t < periods; t++)
ferhat's avatar
ferhat committed
542
    {
543
544
545
546
547
      ti_y_kmin = -min(t, y_kmin);
      ti_y_kmax = min(periods-(t+1), y_kmax);
      it4 = IM.begin();
      eq = -1;
      while (it4 != IM.end())
ferhat's avatar
ferhat committed
548
        {
549
550
551
552
553
          var = it4->first.first.second;
          if (eq != it4->first.first.first+Size*t)
            tmp_b = 0;
          eq = it4->first.first.first+Size*t;
          if (var < (periods+y_kmax)*Size)
ferhat's avatar
ferhat committed
554
            {
555
556
              lag = it4->first.second;
              if (lag <= ti_y_kmax && lag >= ti_y_kmin)
ferhat's avatar
ferhat committed
557
                {
558
                  var += Size*t;
ferhat's avatar
ferhat committed
559
560
561
                }
              else
                {
562
                  tmp_b += u[it4->second+u_count_init*t]*y[index_vara[var+Size*(y_kmin+t)]];
ferhat's avatar
ferhat committed
563
564
565
566
                }
            }
          else
            {
567
568
              b[eq] = it4->second+u_count_init*t;
              u[b[eq]] += tmp_b;
ferhat's avatar
ferhat committed
569
570
571
572
573
574
            }
          it4++;
        }
    }
}

575
576
int
SparseMatrix::Get_u()
ferhat's avatar
ferhat committed
577
578
579
{
  if (!u_liste.empty())
    {
580
      int i = u_liste.back();
ferhat's avatar
ferhat committed
581
582
583
584
585
      u_liste.pop_back();
      return i;
    }
  else
    {
586
      if (u_count < u_count_alloc)
ferhat's avatar
ferhat committed
587
        {
588
          int i = u_count;
ferhat's avatar
ferhat committed
589
590
591
592
593
          u_count++;
          return i;
        }
      else
        {
594
595
          u_count_alloc += 5*u_count_alloc_save;
          u = (double *) mxRealloc(u, u_count_alloc*sizeof(double));
ferhat's avatar
ferhat committed
596
597
          if (!u)
            {
598
              mexPrintf("Error in Get_u: memory exhausted (realloc(%d))\n", u_count_alloc*sizeof(double));
ferhat's avatar
ferhat committed
599
600
601
              mexEvalString("st=fclose('all');clear all;");
              mexErrMsgTxt("Exit from Dynare");
            }
602
          int i = u_count;
ferhat's avatar
ferhat committed
603
604
605
606
607
608
          u_count++;
          return i;
        }
    }
}

609
610
void
SparseMatrix::Delete_u(int pos)
ferhat's avatar
ferhat committed
611
612
613
614
{
  u_liste.push_back(pos);
}

615
616
void
SparseMatrix::Clear_u()
ferhat's avatar
ferhat committed
617
618
619
620
{
  u_liste.clear();
}

621
622
void
SparseMatrix::Print_u()
ferhat's avatar
ferhat committed
623
{
624
625
  for (unsigned int i = 0; i < u_liste.size(); i++)
    mexPrintf("%d ", u_liste[i]);
ferhat's avatar
ferhat committed
626
627
}

628
629
void
SparseMatrix::End(int Size)
ferhat's avatar
ferhat committed
630
631
632
633
634
635
636
637
638
{
  mem_mngr.Free_All();
  mxFree(FNZE_R);
  mxFree(FNZE_C);
  mxFree(NbNZRow);
  mxFree(NbNZCol);
  mxFree(b);
  mxFree(line_done);
  mxFree(pivot);
639
  mxFree(pivot_save);
ferhat's avatar
ferhat committed
640
641
642
643
644
645
  mxFree(pivotk);
  mxFree(pivotv);
  mxFree(pivotva);
}

bool
646
SparseMatrix::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
647
{
648
649
650
651
  long int i, j, nop = nop4/2, t, k;
  double r = 0.0;
  bool OK = true;
  clock_t t001;
ferhat's avatar
ferhat committed
652
653
  t_save_op_s *save_op_s, *save_opa_s, *save_opaa_s;
  int *diff1, *diff2;
654
655
656
657
658
659
  t001 = clock();
  diff1 = (int *) mxMalloc(nop*sizeof(int));
  diff2 = (int *) mxMalloc(nop*sizeof(int));
  int max_save_ops_first = -1;
  j = k = i = 0;
  while (i < nop4 && OK)
ferhat's avatar
ferhat committed
660
    {
661
662
663
664
665
      save_op_s = (t_save_op_s *) &(save_op[i]);
      save_opa_s = (t_save_op_s *) &(save_opa[i]);
      save_opaa_s = (t_save_op_s *) &(save_opaa[i]);
      diff1[j] = save_op_s->first-save_opa_s->first;
      if (max_save_ops_first < save_op_s->first+diff1[j]*(periods-beg_t))
666
        {
667
          max_save_ops_first = save_op_s->first+diff1[j]*(periods-beg_t);
668
        }
ferhat's avatar
ferhat committed
669
670
      switch (save_op_s->operat)
        {
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
        case IFLD:
        case IFDIV:
          OK = (save_op_s->operat == save_opa_s->operat && save_opa_s->operat == save_opaa_s->operat
                && diff1[j] == (save_opa_s->first-save_opaa_s->first));
          i += 2;
          break;
        case IFLESS:
        case IFSUB:
          diff2[j] = save_op_s->second-save_opa_s->second;
          OK = (save_op_s->operat == save_opa_s->operat && save_opa_s->operat == save_opaa_s->operat
                && diff1[j] == (save_opa_s->first-save_opaa_s->first)
                && diff2[j] == (save_opa_s->second-save_opaa_s->second));
          i += 3;
          break;
        default:
          mexPrintf("unknown operator = %d ", save_op_s->operat);
          mexEvalString("st=fclose('all');clear all;");
          filename += " stopped";
          mexErrMsgTxt(filename.c_str());
          break;
ferhat's avatar
ferhat committed
691
692
693
694
695
696
        }
      j++;
    }
  // the same pivot for all remaining periods
  if (OK)
    {
697
698
699
700
701
702
703
704
705
706
      //#pragma omp parallel for  num_threads(atoi(getenv("DYNARE_NUM_THREADS"))) ordered private(j) schedule(dynamic)
      for (i = beg_t; i < periods; i++)
        {
          for (j = 0; j < Size; j++)
            {
              ///#pragma omp ordered
              pivot[i*Size+j] = pivot[(i-1)*Size+j]+Size;
            }
        }
      if (max_save_ops_first >= u_count_alloc)
ferhat's avatar
ferhat committed
707
        {
708
709
          u_count_alloc += max_save_ops_first;
          u = (double *) mxRealloc(u, u_count_alloc*sizeof(double));
710
711
          if (!u)
            {
712
              mexPrintf("Error in Get_u: memory exhausted (realloc(%d))\n", u_count_alloc*sizeof(double));
713
714
715
              mexEvalString("st=fclose('all');clear all;");
              mexErrMsgTxt("Exit from Dynare");
            }
716
717
718
        }
      double *up;
      for (t = 1; t < periods-beg_t-y_kmax; t++)
719
        {
720
721
          i = j = 0;
          while (i < nop4)
ferhat's avatar
ferhat committed
722
            {
723
724
              save_op_s = (t_save_op_s *) (&(save_op[i]));
              up = &u[save_op_s->first+t*diff1[j]];
ferhat's avatar
ferhat committed
725
726
              switch (save_op_s->operat)
                {
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
                case IFLD:
                  r = *up;
                  i += 2;
                  break;
                case IFDIV:
                  *up /= r;
                  i += 2;
                  break;
                case IFSUB:
                  *up -= u[save_op_s->second+t*diff2[j]]*r;;
                  i += 3;
                  break;
                case IFLESS:
                  *up = -u[save_op_s->second+t*diff2[j]]*r;
                  i += 3;
                  break;
ferhat's avatar
ferhat committed
743
744
                }
              j++;
745
            }
ferhat's avatar
ferhat committed
746
        }
747
748
749
      int t1 = max(1, periods-beg_t-y_kmax);
      int periods_beg_t = periods-beg_t;
      for (t = t1; t < periods_beg_t; t++)
ferhat's avatar
ferhat committed
750
        {
751
752
          i = j = 0;
          while (i < nop4)
ferhat's avatar
ferhat committed
753
            {
754
755
              save_op_s = (t_save_op_s *) (&(save_op[i]));
              if (save_op_s->lag < (periods_beg_t-t))
ferhat's avatar
ferhat committed
756
                {
757
                  up = &u[save_op_s->first+t*diff1[j]];
ferhat's avatar
ferhat committed
758
759
                  switch (save_op_s->operat)
                    {
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
                    case IFLD:
                      r = *up;
                      i += 2;
                      break;
                    case IFDIV:
                      *up /= r;
                      i += 2;
                      break;
                    case IFSUB:
                      *up -= u[save_op_s->second+t*diff2[j]]*r;
                      i += 3;
                      break;
                    case IFLESS:
                      *up = -u[save_op_s->second+t*diff2[j]]*r;
                      i += 3;
                      break;
ferhat's avatar
ferhat committed
776
777
778
779
780
781
                    }
                }
              else
                {
                  switch (save_op_s->operat)
                    {
782
783
784
785
786
787
788
789
                    case IFLD:
                    case IFDIV:
                      i += 2;
                      break;
                    case IFSUB:
                    case IFLESS:
                      i += 3;
                      break;
ferhat's avatar
ferhat committed
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
                    }
                }
              j++;
            }
        }
    }
  mxFree(diff1);
  mxFree(diff2);
  return OK;
}

int
SparseMatrix::complete(int beg_t, int Size, int periods, int *b)
{
  long int i, j, k, nop, nopa, nop1, cal_y, nb_var, pos, t, ti, max_var, min_var;
  NonZeroElem *first;
  int *save_code;
  int *diff;
808
  double yy = 0.0, err;
ferhat's avatar
ferhat committed
809

810
811
812
813
814
  int size_of_save_code = (1+y_kmax)*Size*(Size+1+4)/2*4;
  save_code = (int *) mxMalloc(size_of_save_code*sizeof(int));
  int size_of_diff = (1+y_kmax)*Size*(Size+1+4);
  diff = (int *) mxMalloc(size_of_diff*sizeof(int));
  cal_y = y_size*y_kmin;
ferhat's avatar
ferhat committed
815

816
817
818
  i = (beg_t+1)*Size-1;
  nop = 0;
  for (j = i; j > i-Size; j--)
ferhat's avatar
ferhat committed
819
    {
820
821
822
      pos = pivot[j];
      nb_var = At_Row(pos, &first);
      first = first->NZE_R_N;
ferhat's avatar
ferhat committed
823
      nb_var--;
824
825
826
827
828
829
830
831
      save_code[nop] = IFLDZ;
      save_code[nop+1] = 0;
      save_code[nop+2] = 0;
      save_code[nop+3] = 0;
      if ((nop+3) >= size_of_save_code)
        mexPrintf("out of save_code[%d] (bound=%d)\n", nop+2, size_of_save_code);
      nop += 4;
      for (k = 0; k < nb_var; k++)
ferhat's avatar
ferhat committed
832
        {
833
834
835
836
837
838
839
840
          save_code[nop] = IFMUL;
          save_code[nop+1] = index_vara[first->c_index]+cal_y;
          save_code[nop+2] = first->u_index;
          save_code[nop+3] = first->lag_index;
          if ((nop+3) >= size_of_save_code)
            mexPrintf("out of save_code[%d] (bound=%d)\n", nop+2, size_of_save_code);
          nop += 4;
          first = first->NZE_R_N;
ferhat's avatar
ferhat committed
841
        }
842
843
844
845
846
847
848
849
850
851
852
853
854
855
      save_code[nop] = IFADD;
      save_code[nop+1] = b[pos];
      save_code[nop+2] = 0;
      save_code[nop+3] = 0;
      if ((nop+3) >= size_of_save_code)
        mexPrintf("out of save_code[%d] (bound=%d)\n", nop+2, size_of_save_code);
      nop += 4;
      save_code[nop] = IFSTP;
      save_code[nop+1] = index_vara[j]+y_size*y_kmin;
      save_code[nop+2] = 0;
      save_code[nop+3] = 0;
      if ((nop+2) >= size_of_save_code)
        mexPrintf("out of save_code[%d] (bound=%d)\n", nop+2, size_of_save_code);
      nop += 4;
ferhat's avatar
ferhat committed
856
    }
857
858
859
  i = beg_t*Size-1;
  nop1 = nopa = 0;
  for (j = i; j > i-Size; j--)
ferhat's avatar
ferhat committed
860
    {
861
862
863
      pos = pivot[j];
      nb_var = At_Row(pos, &first);
      first = first->NZE_R_N;
ferhat's avatar
ferhat committed
864
      nb_var--;
865
866
867
868
869
      diff[nopa] = 0;
      diff[nopa+1] = 0;
      nopa += 2;
      nop1 += 4;
      for (k = 0; k < nb_var; k++)
ferhat's avatar
ferhat committed
870
        {
871
872
873
874
875
876
877
878
879
          diff[nopa] = save_code[nop1+1]-(index_vara[first->c_index]+cal_y);
          diff[nopa+1] = save_code[nop1+2]-(first->u_index);
          if ((nop1+2) >= size_of_save_code)
            mexPrintf("out of save_code[%d] (bound=%d)\n", nop1+2, size_of_save_code);
          if ((nopa+1) >= size_of_diff)
            mexPrintf("out of diff[%d] (bound=%d)\n", nopa+2, size_of_diff);
          nopa += 2;
          nop1 += 4;
          first = first->NZE_R_N;
ferhat's avatar
ferhat committed
880
        }
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
      diff[nopa] = save_code[nop1+1]-(b[pos]);
      diff[nopa+1] = 0;
      if ((nop1+3) >= size_of_save_code)
        mexPrintf("out of save_code[%d] (bound=%d)\n", nop1+2, size_of_save_code);
      if ((nopa+1) >= size_of_diff)
        mexPrintf("out of diff[%d] (bound=%d)\n", nopa+2, size_of_diff);
      nopa += 2;
      nop1 += 4;
      diff[nopa] = save_code[nop1+1]-(index_vara[j]+y_size*y_kmin);
      diff[nopa+1] = 0;
      if ((nop1+4) >= size_of_save_code)
        mexPrintf("out of save_code[%d] (bound=%d)\n", nop1+2, size_of_save_code);
      if ((nopa+1) >= size_of_diff)
        mexPrintf("out of diff[%d] (bound=%d)\n", nopa+2, size_of_diff);
      nopa += 2;
      nop1 += 4;
ferhat's avatar
ferhat committed
897
    }
898
899
900
901
  max_var = (periods+y_kmin)*y_size;
  min_var = y_kmin*y_size;
  int k1 = 0;
  for (t = periods+y_kmin-1; t >= beg_t+y_kmin; t--)
ferhat's avatar
ferhat committed
902
    {
903
904
905
      j = 0;
      ti = t-y_kmin-beg_t;
      for (i = 0; i < nop; i += 4)
ferhat's avatar
ferhat committed
906
907
908
        {
          switch (save_code[i])
            {
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
            case IFLDZ:
              yy = 0;
              break;
            case IFMUL:
              k = save_code[i+1]+ti*diff[j];
              if (k < max_var && k > min_var)
                {
                  yy += y[k]*u[save_code[i+2]+ti*diff[j+1]];
                }
              break;
            case IFADD:
              yy = -(yy+u[save_code[i+1]+ti*diff[j]]);
              break;
            case IFSTP:
              k = save_code[i+1]+ti*diff[j];
              k1 = k;
              err = yy - y[k];
              y[k] += slowc*(err);
              break;
ferhat's avatar
ferhat committed
928
            }
929
          j += 2;
ferhat's avatar
ferhat committed
930
931
932
933
        }
    }
  mxFree(save_code);
  mxFree(diff);
934
  return (beg_t);
ferhat's avatar
ferhat committed
935
936
}

937
double
938
SparseMatrix::bksub(int tbreak, int last_period, int Size, double slowc_l)
939
940
941
942
943
{
  NonZeroElem *first;
  int i, j, k;
  double yy;
  res1 = res2 = max_res = 0;
944
945
  for (i = 0; i < y_size*(periods+y_kmin); i++)
    y[i] = ya[i];
946
  if (symbolic && tbreak)
947
    last_period = complete(tbreak, Size, periods, b);
948
  else
949
950
    last_period = periods;
  for (int t = last_period+y_kmin-1; t >= y_kmin; t--)
951
    {
952
953
954
955
      int ti = (t-y_kmin)*Size;
      int cal = y_kmin*Size;
      int cal_y = y_size*y_kmin;
      for (i = ti-1; i >= ti-Size; i--)
956
        {
957
958
959
960
          j = i+cal;
          int pos = pivot[i+Size];
          int nb_var = At_Row(pos, &first);
          first = first->NZE_R_N;
961
          nb_var--;
962
963
964
965
966
          int eq = index_vara[j]+y_size;
          yy = 0;
          ostringstream tmp_out;
          tmp_out.clear();
          for (k = 0; k < nb_var; k++)
967
            {
968
969
970
              yy += y[index_vara[first->c_index]+cal_y]*u[first->u_index];
              tmp_out << "y[" << index_vara[first->c_index]+cal_y << "]" << "(" << y[index_vara[first->c_index]+cal_y] << ")*u[" << first->u_index << "](" << u[first->u_index] << ")+";
              first = first->NZE_R_N;
971
            }
972
973
          yy = -(yy+y[eq]+u[b[pos]]);
          direction[eq] = yy;
974
975
976
977
978
979
980
981
982
          y[eq] += slowc_l*yy;
        }
    }
  return res1;
}

double
SparseMatrix::simple_bksub(int it_, int Size, double slowc_l)
{
983
  int i, k;
984
985
986
  double yy;
  NonZeroElem *first;
  res1 = res2 = max_res = 0;
987
988
989
  for (i = 0; i < y_size; i++)
    y[i+it_*y_size] = ya[i+it_*y_size];
  for (i = Size-1; i >= 0; i--)
990
    {
991
992
993
      int pos = pivot[i];
      int nb_var = At_Row(pos, &first);
      first = first->NZE_R_N;
994
      nb_var--;
995
      int eq = index_vara[i];
996
      yy = 0;
997
998
999
1000
      ostringstream tmp_out;
      tmp_out.clear();
      tmp_out << "|";
      for (k = 0; k < nb_var; k++)