2008-01-11 14:42:14 +01:00
/*
2012-06-08 18:22:34 +02:00
* Copyright ( C ) 2007 - 2012 Dynare Team
2008-01-11 14:42:14 +01:00
*
* 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/>.
*/
2010-02-05 18:50:57 +01:00
//#define _GLIBCXX_USE_C99_FENV_TR1 1
//#include <cfenv>
2010-02-05 12:05:21 +01:00
2008-06-28 13:20:45 +02:00
# include <cstring>
2009-08-25 11:43:01 +02:00
# include <ctime>
2009-07-21 17:50:12 +02:00
# include <sstream>
2007-11-21 00:27:30 +01:00
# include "SparseMatrix.hh"
SparseMatrix : : SparseMatrix ( )
{
2009-08-25 11:43:01 +02:00
pivotva = NULL ;
g_save_op = NULL ;
g_nop_all = 0 ;
2007-11-21 00:27:30 +01:00
mem_mngr . init_Mem ( ) ;
2009-08-25 11:43:01 +02:00
symbolic = true ;
alt_symbolic = false ;
alt_symbolic_count = 0 ;
max_u = 0 ;
min_u = 0x7FFFFFFF ;
res1a = 9.0e60 ;
tbreak_g = 0 ;
start_compare = 0 ;
2009-07-21 17:50:12 +02:00
restart = 0 ;
2009-10-30 17:29:16 +01:00
IM_i . clear ( ) ;
2010-07-23 11:20:24 +02:00
lu_inc_tol = 1e-10 ;
2007-11-21 00:27:30 +01:00
}
2009-08-25 11:43:01 +02:00
int
SparseMatrix : : NRow ( int r )
2007-11-21 00:27:30 +01:00
{
return NbNZRow [ r ] ;
}
2009-08-25 11:43:01 +02:00
int
SparseMatrix : : NCol ( int c )
2007-11-21 00:27:30 +01:00
{
return NbNZCol [ c ] ;
}
2009-08-25 11:43:01 +02:00
int
SparseMatrix : : At_Row ( int r , NonZeroElem * * first )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
( * first ) = FNZE_R [ r ] ;
2007-11-21 00:27:30 +01:00
return NbNZRow [ r ] ;
}
int
SparseMatrix : : Union_Row ( int row1 , int row2 )
{
NonZeroElem * first1 , * first2 ;
2009-08-25 11:43:01 +02:00
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 )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
if ( first1 - > c_index = = first2 - > c_index )
2007-11-21 00:27:30 +01:00
{
nb_elem + + ;
i1 + + ;
i2 + + ;
2009-08-25 11:43:01 +02:00
first1 = first1 - > NZE_R_N ;
first2 = first2 - > NZE_R_N ;
2007-11-21 00:27:30 +01:00
}
2009-08-25 11:43:01 +02:00
else if ( first1 - > c_index < first2 - > c_index )
2007-11-21 00:27:30 +01:00
{
nb_elem + + ;
i1 + + ;
2009-08-25 11:43:01 +02:00
first1 = first1 - > NZE_R_N ;
2007-11-21 00:27:30 +01:00
}
else
{
nb_elem + + ;
i2 + + ;
2009-08-25 11:43:01 +02:00
first2 = first2 - > NZE_R_N ;
2007-11-21 00:27:30 +01:00
}
}
return nb_elem ;
}
int
SparseMatrix : : At_Pos ( int r , int c , NonZeroElem * * first )
{
2009-08-25 11:43:01 +02:00
( * first ) = FNZE_R [ r ] ;
while ( ( * first ) - > c_index ! = c )
( * first ) = ( * first ) - > NZE_R_N ;
2007-11-21 00:27:30 +01:00
return NbNZRow [ r ] ;
}
2009-08-25 11:43:01 +02:00
int
SparseMatrix : : At_Col ( int c , NonZeroElem * * first )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
( * first ) = FNZE_C [ c ] ;
2007-11-21 00:27:30 +01:00
return NbNZCol [ c ] ;
}
2009-08-25 11:43:01 +02:00
int
SparseMatrix : : At_Col ( int c , int lag , NonZeroElem * * first )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
( * first ) = FNZE_C [ c ] ;
int i = 0 ;
while ( ( * first ) - > lag_index ! = lag & & ( * first ) )
( * first ) = ( * first ) - > NZE_C_N ;
2007-11-21 00:27:30 +01:00
if ( ( * first ) )
{
2009-08-25 11:43:01 +02:00
NonZeroElem * firsta = ( * first ) ;
2007-11-21 00:27:30 +01:00
if ( ! firsta - > NZE_C_N )
i + + ;
else
{
2009-08-25 11:43:01 +02:00
while ( firsta - > lag_index = = lag & & firsta - > NZE_C_N )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
firsta = firsta - > NZE_C_N ;
2007-11-21 00:27:30 +01:00
i + + ;
}
2011-02-04 16:53:12 +01:00
if ( firsta - > lag_index = = lag )
i + + ;
2007-11-21 00:27:30 +01:00
}
}
return i ;
}
2009-08-25 11:43:01 +02:00
void
SparseMatrix : : Delete ( const int r , const int c )
2007-11-21 00:27:30 +01:00
{
2009-12-16 18:18:38 +01:00
NonZeroElem * first = FNZE_R [ r ] , * firsta = NULL ;
2009-08-25 11:43:01 +02:00
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 ;
2007-11-21 00:27:30 +01:00
NbNZRow [ r ] - - ;
2009-08-25 11:43:01 +02:00
first = FNZE_C [ c ] ;
firsta = NULL ;
while ( first - > r_index ! = r )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
firsta = first ;
first = first - > NZE_C_N ;
2007-11-21 00:27:30 +01:00
}
2009-08-25 11:43:01 +02:00
if ( firsta ! = NULL )
firsta - > NZE_C_N = first - > NZE_C_N ;
if ( first = = FNZE_C [ c ] )
2009-12-16 18:18:38 +01:00
FNZE_C [ c ] = first - > NZE_C_N ;
2009-08-25 11:43:01 +02:00
2007-11-21 00:27:30 +01:00
u_liste . push_back ( first - > u_index ) ;
mem_mngr . mxFree_NZE ( first ) ;
NbNZCol [ c ] - - ;
}
2009-08-25 11:43:01 +02:00
void
SparseMatrix : : Print ( int Size , int * b )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
int a , i , j , k , l ;
2007-11-21 00:27:30 +01:00
mexPrintf ( " " ) ;
2009-08-25 11:43:01 +02:00
for ( k = 0 ; k < Size * periods ; k + + )
mexPrintf ( " %-2d " , k ) ;
2007-11-21 00:27:30 +01:00
mexPrintf ( " | " ) ;
2009-08-25 11:43:01 +02:00
for ( k = 0 ; k < Size * periods ; k + + )
mexPrintf ( " %8d " , k ) ;
2007-11-21 00:27:30 +01:00
mexPrintf ( " \n " ) ;
2009-08-25 11:43:01 +02:00
for ( i = 0 ; i < Size * periods ; i + + )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
NonZeroElem * first = FNZE_R [ i ] ;
j = NbNZRow [ i ] ;
mexPrintf ( " %-2d " , i ) ;
a = 0 ;
for ( k = 0 ; k < j ; k + + )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
for ( l = 0 ; l < ( first - > c_index - a ) ; l + + )
2007-11-21 00:27:30 +01:00
mexPrintf ( " " ) ;
2009-08-25 11:43:01 +02:00
mexPrintf ( " %-2d " , first - > u_index ) ;
a = first - > c_index + 1 ;
first = first - > NZE_R_N ;
2007-11-21 00:27:30 +01:00
}
2009-08-25 11:43:01 +02:00
for ( k = a ; k < Size * periods ; k + + )
2007-11-21 00:27:30 +01:00
mexPrintf ( " " ) ;
2009-08-25 11:43:01 +02:00
mexPrintf ( " %-2d " , b [ i ] ) ;
2007-11-21 00:27:30 +01:00
2009-08-25 11:43:01 +02:00
first = FNZE_R [ i ] ;
j = NbNZRow [ i ] ;
mexPrintf ( " | %-2d " , i ) ;
a = 0 ;
for ( k = 0 ; k < j ; k + + )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
for ( l = 0 ; l < ( first - > c_index - a ) ; l + + )
2007-11-21 00:27:30 +01:00
mexPrintf ( " " ) ;
2009-08-25 11:43:01 +02:00
mexPrintf ( " %8.4f " , double ( u [ first - > u_index ] ) ) ;
a = first - > c_index + 1 ;
first = first - > NZE_R_N ;
2007-11-21 00:27:30 +01:00
}
2009-08-25 11:43:01 +02:00
for ( k = a ; k < Size * periods ; k + + )
2007-11-21 00:27:30 +01:00
mexPrintf ( " " ) ;
2009-08-25 11:43:01 +02:00
mexPrintf ( " %8.4f " , double ( u [ b [ i ] ] ) ) ;
2007-11-21 00:27:30 +01:00
mexPrintf ( " \n " ) ;
}
}
2009-08-25 11:43:01 +02:00
void
SparseMatrix : : Insert ( const int r , const int c , const int u_index , const int lag_index )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
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 ) )
2009-01-30 12:36:15 +01:00
{
2009-08-25 11:43:01 +02:00
firsta = first ;
first = a ;
2007-11-21 00:27:30 +01:00
}
2009-08-25 11:43:01 +02:00
firstn - > u_index = u_index ;
firstn - > r_index = r ;
firstn - > c_index = c ;
firstn - > lag_index = lag_index ;
if ( first - > c_index > c )
2009-01-30 12:36:15 +01:00
{
2009-08-25 11:43:01 +02:00
if ( first = = FNZE_R [ r ] )
FNZE_R [ r ] = firstn ;
if ( firsta ! = NULL )
firsta - > NZE_R_N = firstn ;
firstn - > NZE_R_N = first ;
2007-11-21 00:27:30 +01:00
}
2009-12-16 14:21:31 +01:00
else
2009-01-30 12:36:15 +01:00
{
2009-08-25 11:43:01 +02:00
first - > NZE_R_N = firstn ;
firstn - > NZE_R_N = NULL ;
2009-01-30 12:36:15 +01:00
}
NbNZRow [ r ] + + ;
2009-08-25 11:43:01 +02:00
first = FNZE_C [ c ] ;
firsta = NULL ;
while ( first - > r_index < r & & ( a = first - > NZE_C_N ) )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
firsta = first ;
first = a ;
2007-11-21 00:27:30 +01:00
}
2009-08-25 11:43:01 +02:00
if ( first - > r_index > r )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
if ( first = = FNZE_C [ c ] )
FNZE_C [ c ] = firstn ;
if ( firsta ! = NULL )
firsta - > NZE_C_N = firstn ;
firstn - > NZE_C_N = first ;
2007-11-21 00:27:30 +01:00
}
2009-12-16 14:21:31 +01:00
else
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
first - > NZE_C_N = firstn ;
firstn - > NZE_C_N = NULL ;
2007-11-21 00:27:30 +01:00
}
2009-08-25 11:43:01 +02:00
2007-11-21 00:27:30 +01:00
NbNZCol [ c ] + + ;
}
2009-08-25 11:43:01 +02:00
void
2010-07-23 11:20:24 +02:00
SparseMatrix : : 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 )
2007-11-21 00:27:30 +01:00
{
2009-12-16 14:21:31 +01:00
unsigned int eq , var ;
int i , j , lag ;
2009-08-25 11:43:01 +02:00
filename = file_name ;
2007-11-21 00:27:30 +01:00
mem_mngr . fixe_file_name ( file_name ) ;
if ( ! SaveCode . is_open ( ) )
{
2009-12-16 18:18:38 +01:00
if ( steady_state )
2009-08-25 11:43:01 +02:00
SaveCode . open ( ( file_name + " _static.bin " ) . c_str ( ) , ios : : in | ios : : binary ) ;
2009-12-16 18:18:38 +01:00
else
SaveCode . open ( ( file_name + " _dynamic.bin " ) . c_str ( ) , ios : : in | ios : : binary ) ;
2007-11-21 00:27:30 +01:00
if ( ! SaveCode . is_open ( ) )
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
2009-12-16 18:18:38 +01:00
if ( steady_state )
2010-09-24 12:52:58 +02:00
tmp < < " in Read_SparseMatrix, " < < file_name < < " _static.bin cannot be opened \n " ;
2009-12-16 18:18:38 +01:00
else
2010-09-24 12:52:58 +02:00
tmp < < " in Read_SparseMatrix, " < < file_name < < " _dynamic.bin cannot be opened \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2007-11-21 00:27:30 +01:00
}
}
IM_i . clear ( ) ;
2009-12-16 18:18:38 +01:00
if ( two_boundaries )
{
2010-07-23 11:20:24 +02:00
if ( stack_solve_algo = = 5 )
2009-08-29 18:04:06 +02:00
{
2010-07-23 11:20:24 +02:00
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 + + )
IM_i [ make_pair ( make_pair ( j , Size * ( periods + y_kmax ) ) , 0 ) ] = j ;
2009-08-29 18:04:06 +02:00
}
2010-10-11 19:21:32 +02:00
else if ( stack_solve_algo > = 0 | | stack_solve_algo < = 4 )
2010-07-23 11:20:24 +02:00
{
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 ( var - lag * Size , - lag ) , eq ) ] = j ;
}
for ( j = 0 ; j < Size ; j + + )
IM_i [ make_pair ( make_pair ( Size * ( periods + y_kmax ) , 0 ) , j ) ] = j ;
}
2009-12-16 18:18:38 +01:00
}
else
{
2010-10-22 16:49:47 +02:00
if ( ( stack_solve_algo = = 5 & & ! steady_state ) | | ( solve_algo = = 5 & & steady_state ) )
2010-07-23 11:20:24 +02:00
{
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 ;
}
}
2011-02-04 16:53:12 +01:00
else if ( ( ( stack_solve_algo > = 0 | | stack_solve_algo < = 4 ) & & ! steady_state ) | | ( ( solve_algo > = 6 | | solve_algo < = 8 ) & & steady_state ) )
2009-08-29 18:04:06 +02:00
{
2010-07-23 11:20:24 +02:00
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 ( var - lag * Size , - lag ) , eq ) ] = j ;
}
2009-08-29 18:04:06 +02:00
}
2009-12-16 18:18:38 +01:00
}
2009-08-25 11:43:01 +02:00
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 )
2010-01-22 11:03:29 +01:00
for ( i = 1 ; i < periods + y_kmin + y_kmax ; i + + )
for ( j = 0 ; j < Size ; j + + )
index_vara [ j + Size * i ] = index_vara [ j + Size * ( i - 1 ) ] + y_size ;
2009-08-25 11:43:01 +02:00
index_equa = ( int * ) mxMalloc ( Size * sizeof ( int ) ) ;
for ( j = 0 ; j < Size ; j + + )
2010-01-22 11:03:29 +01:00
SaveCode . read ( reinterpret_cast < char * > ( & index_equa [ j ] ) , sizeof ( * index_equa ) ) ;
2007-11-21 00:27:30 +01:00
}
2009-08-25 11:43:01 +02:00
void
2012-03-06 11:03:13 +01:00
SparseMatrix : : Simple_Init ( int Size , map < pair < pair < int , int > , int > , int > & IM , bool & zero_solution )
2009-01-20 23:04:37 +01:00
{
int i , eq , var , lag ;
2009-08-25 11:43:01 +02:00
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 ) ) ;
2009-01-30 12:36:15 +01:00
2009-01-20 23:04:37 +01:00
mem_mngr . init_CHUNK_BLCK_SIZE ( u_count ) ;
2009-08-25 11:43:01 +02:00
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 ) ;
it4 = IM . begin ( ) ;
eq = - 1 ;
//#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
for ( i = 0 ; i < Size ; i + + )
2009-01-30 12:36:15 +01:00
{
2009-08-25 11:43:01 +02:00
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 ;
2009-01-30 12:36:15 +01:00
}
2009-08-25 11:43:01 +02:00
int u_count1 = Size ;
while ( it4 ! = IM . end ( ) )
2009-01-20 23:04:37 +01:00
{
2009-08-25 11:43:01 +02:00
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*/
2009-01-20 23:04:37 +01:00
{
NbNZRow [ eq ] + + ;
NbNZCol [ var ] + + ;
2009-08-25 11:43:01 +02:00
first = mem_mngr . mxMalloc_NZE ( ) ;
first - > NZE_C_N = NULL ;
first - > NZE_R_N = NULL ;
2009-12-16 14:21:31 +01:00
first - > u_index = u_count1 ;
2009-08-25 11:43:01 +02:00
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 ;
2009-01-20 23:04:37 +01:00
u_count1 + + ;
}
it4 + + ;
}
2009-08-25 11:43:01 +02:00
//#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
2010-10-18 17:28:21 +02:00
double cum_abs_sum = 0 ;
2009-08-25 11:43:01 +02:00
for ( i = 0 ; i < Size ; i + + )
2010-10-18 17:28:21 +02:00
{
b [ i ] = i ;
cum_abs_sum + = fabs ( u [ i ] ) ;
}
if ( cum_abs_sum < 1e-20 )
zero_solution = true ;
else
zero_solution = false ;
2009-01-20 23:04:37 +01:00
mxFree ( temp_NZE_R ) ;
mxFree ( temp_NZE_C ) ;
2009-12-16 14:21:31 +01:00
u_count = u_count1 ;
2009-01-20 23:04:37 +01:00
}
2009-08-25 11:43:01 +02:00
void
2011-01-03 14:15:12 +01:00
SparseMatrix : : Init_Matlab_Sparse_Simple ( int Size , map < pair < pair < int , int > , int > , int > & IM , mxArray * A_m , mxArray * b_m , bool & zero_solution , mxArray * x0_m )
2010-07-23 11:20:24 +02:00
{
int i , eq , var ;
double * b = mxGetPr ( b_m ) ;
if ( ! b )
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Init_Matlab_Sparse_Simple, can't retrieve b vector \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2010-07-23 11:20:24 +02:00
}
2011-01-03 14:15:12 +01:00
double * x0 = mxGetPr ( x0_m ) ;
if ( ! x0 )
{
ostringstream tmp ;
tmp < < " in Init_Matlab_Sparse_Simple, can't retrieve x0 vector \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
}
2010-07-23 11:20:24 +02:00
mwIndex * Ai = mxGetIr ( A_m ) ;
if ( ! Ai )
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Init_Matlab_Sparse_Simple, can't allocate Ai index vector \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2010-07-23 11:20:24 +02:00
}
mwIndex * Aj = mxGetJc ( A_m ) ;
if ( ! Aj )
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Init_Matlab_Sparse_Simple, can't allocate Aj index vector \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2010-07-23 11:20:24 +02:00
}
double * A = mxGetPr ( A_m ) ;
if ( ! A )
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Init_Matlab_Sparse_Simple, can't retrieve A matrix \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2010-07-23 11:20:24 +02:00
}
map < pair < pair < int , int > , int > , int > : : iterator it4 ;
for ( i = 0 ; i < y_size * ( periods + y_kmin ) ; i + + )
ya [ i ] = y [ i ] ;
2010-12-31 16:29:39 +01:00
# ifdef DEBUG
2010-07-23 11:20:24 +02:00
unsigned int max_nze = mxGetNzmax ( A_m ) ;
2010-09-17 09:57:38 +02:00
# endif
2010-07-23 11:20:24 +02:00
unsigned int NZE = 0 ;
int last_var = 0 ;
2010-10-18 17:28:21 +02:00
double cum_abs_sum = 0 ;
2010-07-23 11:20:24 +02:00
for ( i = 0 ; i < Size ; i + + )
2010-10-18 17:28:21 +02:00
{
b [ i ] = u [ i ] ;
cum_abs_sum + = fabs ( b [ i ] ) ;
2011-01-03 14:15:12 +01:00
x0 [ i ] = y [ i ] ;
2010-10-18 17:28:21 +02:00
}
if ( cum_abs_sum < 1e-20 )
zero_solution = true ;
else
zero_solution = false ;
2010-07-23 11:20:24 +02:00
Aj [ 0 ] = 0 ;
2010-10-18 17:28:21 +02:00
last_var = - 1 ;
2010-07-23 11:20:24 +02:00
it4 = IM . begin ( ) ;
while ( it4 ! = IM . end ( ) )
{
var = it4 - > first . first . first ;
if ( var ! = last_var )
{
Aj [ 1 + last_var ] = NZE ;
last_var = var ;
}
eq = it4 - > first . second ;
int index = it4 - > second ;
2010-12-31 16:29:39 +01:00
# ifdef DEBUG
2011-02-04 16:53:12 +01:00
if ( index < 0 | | index > = u_count_alloc | | index > Size + Size * Size )
2010-07-23 11:20:24 +02:00
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Init_Matlab_Sparse_Simple, index ( " < < index < < " ) out of range for u vector max = " < < Size + Size * Size < < " allocated = " < < u_count_alloc < < " \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2010-07-23 11:20:24 +02:00
}
if ( NZE > = max_nze )
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Init_Matlab_Sparse_Simple, exceeds the capacity of A_m sparse matrix \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2010-07-23 11:20:24 +02:00
}
# endif
A [ NZE ] = u [ index ] ;
Ai [ NZE ] = eq ;
NZE + + ;
2010-12-31 16:29:39 +01:00
# ifdef DEBUG
2010-08-18 13:51:57 +02:00
if ( eq < 0 | | eq > = Size )
2010-07-23 11:20:24 +02:00
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Init_Matlab_Sparse_Simple, index ( " < < eq < < " ) out of range for b vector \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2011-02-04 16:53:12 +01:00
}
2010-08-18 13:51:57 +02:00
if ( var < 0 | | var > = Size )
2010-07-23 11:20:24 +02:00
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Init_Matlab_Sparse_Simple, index ( " < < var < < " ) out of range for index_vara vector \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2010-07-23 11:20:24 +02:00
}
2010-08-18 13:51:57 +02:00
if ( index_vara [ var ] < 0 | | index_vara [ var ] > = y_size )
2010-07-23 11:20:24 +02:00
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Init_Matlab_Sparse_Simple, index ( " < < index_vara [ var ] < < " ) out of range for y vector max= " < < y_size < < " (0) \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2010-07-23 11:20:24 +02:00
}
# endif
it4 + + ;
2011-02-04 16:53:12 +01:00
}
2010-07-23 11:20:24 +02:00
Aj [ Size ] = NZE ;
}
void
2011-01-03 14:15:12 +01:00
SparseMatrix : : 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 , mxArray * x0_m )
2010-07-23 11:20:24 +02:00
{
int t , i , eq , var , lag , ti_y_kmin , ti_y_kmax ;
double * b = mxGetPr ( b_m ) ;
if ( ! b )
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Init_Matlab_Sparse, can't retrieve b vector \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2010-07-23 11:20:24 +02:00
}
2011-01-03 14:15:12 +01:00
double * x0 = mxGetPr ( x0_m ) ;
if ( ! x0 )
{
ostringstream tmp ;
tmp < < " in Init_Matlab_Sparse_Simple, can't retrieve x0 vector \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
}
2010-07-23 11:20:24 +02:00
mwIndex * Ai = mxGetIr ( A_m ) ;
if ( ! Ai )
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Init_Matlab_Sparse, can't allocate Ai index vector \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2010-07-23 11:20:24 +02:00
}
mwIndex * Aj = mxGetJc ( A_m ) ;
if ( ! Aj )
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Init_Matlab_Sparse, can't allocate Aj index vector \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2010-07-23 11:20:24 +02:00
}
double * A = mxGetPr ( A_m ) ;
if ( ! A )
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Init_Matlab_Sparse, can't retrieve A matrix \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2010-07-23 11:20:24 +02:00
}
map < pair < pair < int , int > , int > , int > : : iterator it4 ;
for ( i = 0 ; i < y_size * ( periods + y_kmin ) ; i + + )
ya [ i ] = y [ i ] ;
2010-12-31 16:29:39 +01:00
# ifdef DEBUG
2010-07-23 11:20:24 +02:00
unsigned int max_nze = mxGetNzmax ( A_m ) ;
# endif
unsigned int NZE = 0 ;
int last_var = 0 ;
for ( i = 0 ; i < periods * Size ; i + + )
2011-01-03 14:15:12 +01:00
{
b [ i ] = 0 ;
x0 [ i ] = y [ index_vara [ Size * y_kmin + i ] ] ;
}
2010-07-23 11:20:24 +02:00
Aj [ 0 ] = 0 ;
for ( t = 0 ; t < periods ; t + + )
{
last_var = 0 ;
it4 = IM . begin ( ) ;
while ( it4 ! = IM . end ( ) )
{
var = it4 - > first . first . first ;
if ( var ! = last_var )
{
Aj [ 1 + last_var + t * Size ] = NZE ;
last_var = var ;
}
eq = it4 - > first . second + Size * t ;
lag = - it4 - > first . first . second ;
int index = it4 - > second + ( t - lag ) * u_count_init ;
if ( var < ( periods + y_kmax ) * Size )
{
ti_y_kmin = - min ( t , y_kmin ) ;
2011-02-04 16:53:12 +01:00
ti_y_kmax = min ( periods - ( t + 1 ) , y_kmax ) ;
2010-07-23 11:20:24 +02:00
int ti_new_y_kmax = min ( t , y_kmax ) ;
int ti_new_y_kmin = - min ( periods - ( t + 1 ) , y_kmin ) ;
if ( lag < = ti_new_y_kmax & & lag > = ti_new_y_kmin ) /*Build the index for sparse matrix containing the jacobian : u*/
{
2010-12-31 16:29:39 +01:00
# ifdef DEBUG
2011-02-04 16:53:12 +01:00
if ( index < 0 | | index > = u_count_alloc | | index > Size + Size * Size )
{
ostringstream tmp ;
tmp < < " in Init_Matlab_Sparse, index ( " < < index < < " ) out of range for u vector max = " < < Size + Size * Size < < " allocated = " < < u_count_alloc < < " \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
}
if ( NZE > = max_nze )
{
ostringstream tmp ;
tmp < < " in Init_Matlab_Sparse, exceeds the capacity of A_m sparse matrix \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
}
2010-07-23 11:20:24 +02:00
# endif
A [ NZE ] = u [ index ] ;
Ai [ NZE ] = eq - lag * Size ;
NZE + + ;
}
if ( lag > ti_y_kmax | | lag < ti_y_kmin )
{
2010-12-31 16:29:39 +01:00
# ifdef DEBUG
2010-09-24 12:52:58 +02:00
if ( eq < 0 | | eq > = Size * periods )
2010-07-23 11:20:24 +02:00
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Init_Matlab_Sparse, index ( " < < eq < < " ) out of range for b vector \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2011-02-04 16:53:12 +01:00
}
2010-07-23 11:20:24 +02:00
if ( var + Size * ( y_kmin + t + lag ) < 0 | | var + Size * ( y_kmin + t + lag ) > = Size * ( periods + y_kmin + y_kmax ) )
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Init_Matlab_Sparse, index ( " < < var + Size * ( y_kmin + t + lag ) < < " ) out of range for index_vara vector \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2010-07-23 11:20:24 +02:00
}
if ( index_vara [ var + Size * ( y_kmin + t + lag ) ] < 0 | | index_vara [ var + Size * ( y_kmin + t + lag ) ] > = y_size * ( periods + y_kmin + y_kmax ) )
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Init_Matlab_Sparse, index ( " < < index_vara [ var + Size * ( y_kmin + t + lag ) ] < < " ) out of range for y vector max= " < < y_size * ( periods + y_kmin + y_kmax ) < < " \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2010-07-23 11:20:24 +02:00
}
# endif
b [ eq ] + = u [ index + lag * u_count_init ] * y [ index_vara [ var + Size * ( y_kmin + t + lag ) ] ] ;
}
}
else /* ...and store it in the u vector*/
{
2010-12-31 16:29:39 +01:00
# ifdef DEBUG
2010-07-23 11:20:24 +02:00
if ( index < 0 | | index > = u_count_alloc )
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Init_Matlab_Sparse, index ( " < < index < < " ) out of range for u vector \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2010-07-23 11:20:24 +02:00
}
if ( eq < 0 | | eq > = ( Size * periods ) )
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Init_Matlab_Sparse, index ( " < < eq < < " ) out of range for b vector \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2010-07-23 11:20:24 +02:00
}
# endif
b [ eq ] + = u [ index ] ;
}
it4 + + ;
}
}
Aj [ Size * periods ] = NZE ;
}
void
SparseMatrix : : Init_GE ( int periods , int y_kmin , int y_kmax , int Size , map < pair < pair < int , int > , int > , int > & IM )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
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 ) ) ;
2007-11-21 00:27:30 +01:00
mem_mngr . init_CHUNK_BLCK_SIZE ( u_count ) ;
2009-08-25 11:43:01 +02:00
g_save_op = NULL ;
g_nop_all = 0 ;
2009-12-16 18:18:38 +01:00
i = ( periods + y_kmax + 1 ) * Size * sizeof ( NonZeroElem * ) ;
2009-08-25 11:43:01 +02:00
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 + + )
2009-01-30 12:36:15 +01:00
{
2009-08-25 11:43:01 +02:00
b [ i ] = 0 ;
line_done [ i ] = 0 ;
2009-01-30 12:36:15 +01:00
}
2009-08-25 11:43:01 +02:00
//#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
for ( i = 0 ; i < ( periods + y_kmax + 1 ) * Size ; i + + )
2009-01-30 12:36:15 +01:00
{
2009-08-25 11:43:01 +02:00
FNZE_C [ i ] = 0 ;
FNZE_R [ i ] = 0 ;
temp_NZE_C [ i ] = NULL ;
temp_NZE_R [ i ] = NULL ;
NbNZRow [ i ] = 0 ;
NbNZCol [ i ] = 0 ;
2009-01-30 12:36:15 +01:00
}
2009-08-25 11:43:01 +02:00
//#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 + + )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
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 ( ) )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
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 )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
lag = it4 - > first . second ;
if ( lag < = ti_y_kmax & & lag > = ti_y_kmin ) /*Build the index for sparse matrix containing the jacobian : u*/
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
var + = Size * t ;
2007-11-21 00:27:30 +01:00
NbNZRow [ eq ] + + ;
NbNZCol [ var ] + + ;
2009-08-25 11:43:01 +02:00
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 ;
2007-11-21 00:27:30 +01:00
}
2009-07-21 17:50:12 +02:00
else /*Build the additive terms ooutside the simulation periods related to the first lags and the last leads...*/
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
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 ) ] ] ;
2009-07-21 17:50:12 +02:00
2009-08-25 11:43:01 +02:00
}
2007-11-21 00:27:30 +01:00
}
}
2008-12-19 11:24:31 +01:00
else /* ...and store it in the u vector*/
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
b [ eq ] = it4 - > second + u_count_init * t ;
u [ b [ eq ] ] + = tmp_b ;
2009-07-21 17:50:12 +02:00
tmp_b = 0 ;
2007-11-21 00:27:30 +01:00
}
it4 + + ;
}
}
mxFree ( temp_NZE_R ) ;
mxFree ( temp_NZE_C ) ;
}
2009-08-25 11:43:01 +02:00
int
SparseMatrix : : Get_u ( )
2007-11-21 00:27:30 +01:00
{
if ( ! u_liste . empty ( ) )
{
2009-08-25 11:43:01 +02:00
int i = u_liste . back ( ) ;
2007-11-21 00:27:30 +01:00
u_liste . pop_back ( ) ;
return i ;
}
else
{
2009-08-25 11:43:01 +02:00
if ( u_count < u_count_alloc )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
int i = u_count ;
2007-11-21 00:27:30 +01:00
u_count + + ;
return i ;
}
else
{
2009-08-25 11:43:01 +02:00
u_count_alloc + = 5 * u_count_alloc_save ;
u = ( double * ) mxRealloc ( u , u_count_alloc * sizeof ( double ) ) ;
2007-11-21 00:27:30 +01:00
if ( ! u )
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Get_u, memory exhausted (realloc( " < < u_count_alloc * sizeof ( double ) < < " )) \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2007-11-21 00:27:30 +01:00
}
2009-08-25 11:43:01 +02:00
int i = u_count ;
2007-11-21 00:27:30 +01:00
u_count + + ;
return i ;
}
}
}
2009-08-25 11:43:01 +02:00
void
SparseMatrix : : Delete_u ( int pos )
2007-11-21 00:27:30 +01:00
{
u_liste . push_back ( pos ) ;
}
2009-08-25 11:43:01 +02:00
void
SparseMatrix : : Clear_u ( )
2007-11-21 00:27:30 +01:00
{
u_liste . clear ( ) ;
}
2009-08-25 11:43:01 +02:00
void
SparseMatrix : : Print_u ( )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
for ( unsigned int i = 0 ; i < u_liste . size ( ) ; i + + )
mexPrintf ( " %d " , u_liste [ i ] ) ;
2007-11-21 00:27:30 +01:00
}
2009-08-25 11:43:01 +02:00
void
2010-07-23 11:20:24 +02:00
SparseMatrix : : End_GE ( int Size )
2007-11-21 00:27:30 +01:00
{
mem_mngr . Free_All ( ) ;
mxFree ( FNZE_R ) ;
mxFree ( FNZE_C ) ;
mxFree ( NbNZRow ) ;
mxFree ( NbNZCol ) ;
mxFree ( b ) ;
mxFree ( line_done ) ;
mxFree ( pivot ) ;
2009-07-21 17:50:12 +02:00
mxFree ( pivot_save ) ;
2007-11-21 00:27:30 +01:00
mxFree ( pivotk ) ;
mxFree ( pivotv ) ;
mxFree ( pivotva ) ;
}
bool
2009-08-25 11:43:01 +02:00
SparseMatrix : : compare ( int * save_op , int * save_opa , int * save_opaa , int beg_t , int periods , long int nop4 , int Size )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
long int i , j , nop = nop4 / 2 , t , k ;
double r = 0.0 ;
bool OK = true ;
2007-11-21 00:27:30 +01:00
t_save_op_s * save_op_s , * save_opa_s , * save_opaa_s ;
int * diff1 , * diff2 ;
2009-08-25 11:43:01 +02:00
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 )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
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 ) )
2009-01-30 12:36:15 +01:00
{
2009-08-25 11:43:01 +02:00
max_save_ops_first = save_op_s - > first + diff1 [ j ] * ( periods - beg_t ) ;
2009-01-30 12:36:15 +01:00
}
2007-11-21 00:27:30 +01:00
switch ( save_op_s - > operat )
{
2009-08-25 11:43:01 +02:00
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 :
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in compare, unknown operator = " < < save_op_s - > operat < < " \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2007-11-21 00:27:30 +01:00
}
j + + ;
}
// the same pivot for all remaining periods
if ( OK )
{
2009-08-25 11:43:01 +02:00
//#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 )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
u_count_alloc + = max_save_ops_first ;
u = ( double * ) mxRealloc ( u , u_count_alloc * sizeof ( double ) ) ;
2009-01-30 12:36:15 +01:00
if ( ! u )
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in compare, memory exhausted (realloc( " < < u_count_alloc * sizeof ( double ) < < " )) \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2009-01-30 12:36:15 +01:00
}
2009-08-25 11:43:01 +02:00
}
double * up ;
for ( t = 1 ; t < periods - beg_t - y_kmax ; t + + )
2009-01-30 12:36:15 +01:00
{
2009-08-25 11:43:01 +02:00
i = j = 0 ;
while ( i < nop4 )
2007-11-21 00:27:30 +01:00
{
2011-02-04 16:53:12 +01:00
save_op_s = ( t_save_op_s * ) ( & ( save_op [ i ] ) ) ;
2009-08-25 11:43:01 +02:00
up = & u [ save_op_s - > first + t * diff1 [ j ] ] ;
2007-11-21 00:27:30 +01:00
switch ( save_op_s - > operat )
{
2009-08-25 11:43:01 +02:00
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 ;
2007-11-21 00:27:30 +01:00
}
j + + ;
2009-01-30 12:36:15 +01:00
}
2009-06-05 16:45:23 +02:00
}
2009-08-25 11:43:01 +02:00
int t1 = max ( 1 , periods - beg_t - y_kmax ) ;
int periods_beg_t = periods - beg_t ;
for ( t = t1 ; t < periods_beg_t ; t + + )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
i = j = 0 ;
while ( i < nop4 )
2007-11-21 00:27:30 +01:00
{
2011-02-04 16:53:12 +01:00
save_op_s = ( t_save_op_s * ) ( & ( save_op [ i ] ) ) ;
2009-08-25 11:43:01 +02:00
if ( save_op_s - > lag < ( periods_beg_t - t ) )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
up = & u [ save_op_s - > first + t * diff1 [ j ] ] ;
2007-11-21 00:27:30 +01:00
switch ( save_op_s - > operat )
{
2009-08-25 11:43:01 +02:00
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 ;
2007-11-21 00:27:30 +01:00
}
}
else
{
switch ( save_op_s - > operat )
{
2009-08-25 11:43:01 +02:00
case IFLD :
case IFDIV :
i + = 2 ;
break ;
case IFSUB :
case IFLESS :
i + = 3 ;
break ;
2007-11-21 00:27:30 +01:00
}
}
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 ;
2009-08-25 11:43:01 +02:00
double yy = 0.0 , err ;
2007-11-21 00:27:30 +01:00
2009-08-25 11:43:01 +02:00
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 ;
2007-11-21 00:27:30 +01:00
2009-08-25 11:43:01 +02:00
i = ( beg_t + 1 ) * Size - 1 ;
nop = 0 ;
for ( j = i ; j > i - Size ; j - - )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
pos = pivot [ j ] ;
nb_var = At_Row ( pos , & first ) ;
first = first - > NZE_R_N ;
2007-11-21 00:27:30 +01:00
nb_var - - ;
2009-08-25 11:43:01 +02:00
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 + + )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
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 ;
2007-11-21 00:27:30 +01:00
}
2009-08-25 11:43:01 +02:00
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 ;
2007-11-21 00:27:30 +01:00
}
2009-08-25 11:43:01 +02:00
i = beg_t * Size - 1 ;
nop1 = nopa = 0 ;
for ( j = i ; j > i - Size ; j - - )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
pos = pivot [ j ] ;
nb_var = At_Row ( pos , & first ) ;
first = first - > NZE_R_N ;
2007-11-21 00:27:30 +01:00
nb_var - - ;
2009-08-25 11:43:01 +02:00
diff [ nopa ] = 0 ;
diff [ nopa + 1 ] = 0 ;
nopa + = 2 ;
nop1 + = 4 ;
for ( k = 0 ; k < nb_var ; k + + )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
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 ;
2007-11-21 00:27:30 +01:00
}
2009-08-25 11:43:01 +02:00
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 ;
2007-11-21 00:27:30 +01:00
}
2009-08-25 11:43:01 +02:00
max_var = ( periods + y_kmin ) * y_size ;
min_var = y_kmin * y_size ;
for ( t = periods + y_kmin - 1 ; t > = beg_t + y_kmin ; t - - )
2007-11-21 00:27:30 +01:00
{
2009-08-25 11:43:01 +02:00
j = 0 ;
ti = t - y_kmin - beg_t ;
for ( i = 0 ; i < nop ; i + = 4 )
2007-11-21 00:27:30 +01:00
{
switch ( save_code [ i ] )
{
2009-08-25 11:43:01 +02:00
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 ] ;
err = yy - y [ k ] ;
y [ k ] + = slowc * ( err ) ;
break ;
2007-11-21 00:27:30 +01:00
}
2009-08-25 11:43:01 +02:00
j + = 2 ;
2007-11-21 00:27:30 +01:00
}
}
mxFree ( save_code ) ;
mxFree ( diff ) ;
2009-08-25 11:43:01 +02:00
return ( beg_t ) ;
2007-11-21 00:27:30 +01:00
}
2012-02-17 10:59:39 +01:00
void
2009-08-25 11:43:01 +02:00
SparseMatrix : : bksub ( int tbreak , int last_period , int Size , double slowc_l )
2009-01-20 23:04:37 +01:00
{
NonZeroElem * first ;
int i , j , k ;
double yy ;
2009-08-25 11:43:01 +02:00
for ( i = 0 ; i < y_size * ( periods + y_kmin ) ; i + + )
y [ i ] = ya [ i ] ;
2009-01-20 23:04:37 +01:00
if ( symbolic & & tbreak )
2009-08-25 11:43:01 +02:00
last_period = complete ( tbreak , Size , periods , b ) ;
2009-01-20 23:04:37 +01:00
else
2009-08-25 11:43:01 +02:00
last_period = periods ;
for ( int t = last_period + y_kmin - 1 ; t > = y_kmin ; t - - )
2009-01-20 23:04:37 +01:00
{
2009-08-25 11:43:01 +02:00
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 - - )
2009-01-20 23:04:37 +01:00
{
2009-08-25 11:43:01 +02:00
j = i + cal ;
int pos = pivot [ i + Size ] ;
int nb_var = At_Row ( pos , & first ) ;
first = first - > NZE_R_N ;
2009-01-20 23:04:37 +01:00
nb_var - - ;
2009-08-25 11:43:01 +02:00
int eq = index_vara [ j ] + y_size ;
yy = 0 ;
for ( k = 0 ; k < nb_var ; k + + )
2009-01-20 23:04:37 +01:00
{
2009-08-25 11:43:01 +02:00
yy + = y [ index_vara [ first - > c_index ] + cal_y ] * u [ first - > u_index ] ;
first = first - > NZE_R_N ;
2009-01-20 23:04:37 +01:00
}
2009-08-25 11:43:01 +02:00
yy = - ( yy + y [ eq ] + u [ b [ pos ] ] ) ;
direction [ eq ] = yy ;
2009-01-20 23:04:37 +01:00
y [ eq ] + = slowc_l * yy ;
}
}
}
2012-02-17 10:59:39 +01:00
void
2009-01-20 23:04:37 +01:00
SparseMatrix : : simple_bksub ( int it_ , int Size , double slowc_l )
{
2009-08-25 11:43:01 +02:00
int i , k ;
2009-01-20 23:04:37 +01:00
double yy ;
NonZeroElem * first ;
2009-08-25 11:43:01 +02:00
for ( i = 0 ; i < y_size ; i + + )
y [ i + it_ * y_size ] = ya [ i + it_ * y_size ] ;
for ( i = Size - 1 ; i > = 0 ; i - - )
2009-01-20 23:04:37 +01:00
{
2009-08-25 11:43:01 +02:00
int pos = pivot [ i ] ;
int nb_var = At_Row ( pos , & first ) ;
first = first - > NZE_R_N ;
2009-01-20 23:04:37 +01:00
nb_var - - ;
2009-08-25 11:43:01 +02:00
int eq = index_vara [ i ] ;
2009-01-20 23:04:37 +01:00
yy = 0 ;
2009-08-25 11:43:01 +02:00
for ( k = 0 ; k < nb_var ; k + + )
2009-01-20 23:04:37 +01:00
{
2009-08-25 11:43:01 +02:00
yy + = y [ index_vara [ first - > c_index ] + it_ * y_size ] * u [ first - > u_index ] ;
first = first - > NZE_R_N ;
2009-01-20 23:04:37 +01:00
}
2009-08-25 11:43:01 +02:00
yy = - ( yy + y [ eq + it_ * y_size ] + u [ b [ pos ] ] ) ;
direction [ eq + it_ * y_size ] = yy ;
2009-01-20 23:04:37 +01:00
y [ eq + it_ * y_size ] + = slowc_l * yy ;
}
}
2010-09-24 12:52:58 +02:00
void
SparseMatrix : : CheckIt ( int y_size , int y_kmin , int y_kmax , int Size , int periods , int iter )
2009-01-20 23:04:37 +01:00
{
2010-09-24 12:52:58 +02:00
const double epsilon = 1e-7 ;
fstream SaveResult ;
ostringstream out ;
out < < " Result " < < iter ;
SaveResult . open ( out . str ( ) . c_str ( ) , ios : : in ) ;
if ( ! SaveResult . is_open ( ) )
2009-01-20 23:04:37 +01:00
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in CheckIt, Result file cannot be opened \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
}
mexPrintf ( " Reading Result... " ) ;
int row , col ;
SaveResult > > row ;
mexPrintf ( " row=%d \n " , row ) ;
SaveResult > > col ;
mexPrintf ( " col=%d \n " , col ) ;
double G1a ;
mexPrintf ( " Allocated \n " ) ;
NonZeroElem * first ;
for ( int j = 0 ; j < col ; j + + )
{
mexPrintf ( " j=%d " , j ) ;
int nb_equ = At_Col ( j , & first ) ;
mexPrintf ( " nb_equ=%d \n " , nb_equ ) ;
int line ;
if ( first )
line = first - > r_index ;
else
line = - 9999999 ;
for ( int i = 0 ; i < row ; i + + )
2009-08-25 11:43:01 +02:00
{
2010-09-24 12:52:58 +02:00
SaveResult > > G1a ;
if ( line = = i )
2009-08-29 18:04:06 +02:00
{
2010-09-24 12:52:58 +02:00
if ( abs ( u [ first - > u_index ] / G1a - 1 ) > epsilon )
mexPrintf ( " Problem at r=%d c=%d u[first->u_index]=%5.14f G1a[i][j]=%5.14f %f \n " , i , j , u [ first - > u_index ] , G1a , u [ first - > u_index ] / G1a - 1 ) ;
first = first - > NZE_C_N ;
if ( first )
line = first - > r_index ;
2010-01-22 11:03:29 +01:00
else
2010-09-24 12:52:58 +02:00
line = - 9999999 ;
2010-02-05 18:33:34 +01:00
}
else
{
2010-09-24 12:52:58 +02:00
if ( G1a ! = 0.0 )
mexPrintf ( " Problem at r=%d c=%d G1a[i][j]=%f \n " , i , j , G1a ) ;
2010-02-05 18:33:34 +01:00
}
}
2009-01-20 23:04:37 +01:00
}
2010-09-24 12:52:58 +02:00
mexPrintf ( " G1a red done \n " ) ;
SaveResult > > row ;
mexPrintf ( " row(2)=%d \n " , row ) ;
double * B ;
B = ( double * ) mxMalloc ( row * sizeof ( double ) ) ;
for ( int i = 0 ; i < row ; i + + )
SaveResult > > B [ i ] ;
SaveResult . close ( ) ;
mexPrintf ( " done \n " ) ;
mexPrintf ( " Comparing... " ) ;
for ( int i = 0 ; i < row ; i + + )
2009-12-16 18:18:38 +01:00
{
2010-09-24 12:52:58 +02:00
if ( abs ( u [ b [ i ] ] + B [ i ] ) > epsilon )
mexPrintf ( " Problem at i=%d u[b[i]]=%f B[i]=%f \n " , i , u [ b [ i ] ] , B [ i ] ) ;
2009-08-29 18:04:06 +02:00
}
2010-09-24 12:52:58 +02:00
mxFree ( B ) ;
}
void
SparseMatrix : : Check_the_Solution ( int periods , int y_kmin , int y_kmax , int Size , double * u , int * pivot , int * b )
{
const double epsilon = 1e-10 ;
Init_GE ( periods , y_kmin , y_kmax , Size , IM_i ) ;
NonZeroElem * first ;
int cal_y = y_kmin * Size ;
mexPrintf ( " " ) ;
for ( int i = 0 ; i < Size ; i + + )
mexPrintf ( " %8d " , i ) ;
mexPrintf ( " \n " ) ;
for ( int t = y_kmin ; t < periods + y_kmin ; t + + )
2010-02-05 18:33:34 +01:00
{
2010-09-24 12:52:58 +02:00
mexPrintf ( " t=%5d " , t ) ;
for ( int i = 0 ; i < Size ; i + + )
mexPrintf ( " %d %1.6f " , t * y_size + index_vara [ i ] , y [ t * y_size + index_vara [ i ] ] ) ;
mexPrintf ( " \n " ) ;
2010-02-05 18:33:34 +01:00
}
2010-09-24 12:52:58 +02:00
for ( int i = 0 ; i < Size * periods ; i + + )
2010-07-23 11:20:24 +02:00
{
2010-09-24 12:52:58 +02:00
double res = 0 ;
int pos = pivot [ i ] ;
mexPrintf ( " pos[%d]=%d " , i , pos ) ;
int nb_var = At_Row ( pos , & first ) ;
mexPrintf ( " nb_var=%d \n " , nb_var ) ;
for ( int j = 0 ; j < nb_var ; j + + )
2010-07-23 11:20:24 +02:00
{
2010-09-24 12:52:58 +02:00
mexPrintf ( " (y[%d]=%f)*(u[%d]=%f)(r=%d, c=%d) \n " , index_vara [ first - > c_index ] + cal_y , y [ index_vara [ first - > c_index ] + cal_y ] , first - > u_index , u [ first - > u_index ] , first - > r_index , first - > c_index ) ;
res + = y [ index_vara [ first - > c_index ] + cal_y ] * u [ first - > u_index ] ;
first = first - > NZE_R_N ;
2010-07-23 11:20:24 +02:00
}
2010-09-24 12:52:58 +02:00
double tmp_ = res ;
res + = u [ b [ pos ] ] ;
if ( abs ( res ) > epsilon )
mexPrintf ( " Error for equation %d => res=%f y[%d]=%f u[b[%d]]=%f somme(y*u)=%f \n " , pos , res , pos , y [ index_vara [ pos ] ] , pos , u [ b [ pos ] ] , tmp_ ) ;
2010-07-23 11:20:24 +02:00
}
2010-09-24 12:52:58 +02:00
}
2011-02-04 16:53:12 +01:00
mxArray *
SparseMatrix : : substract_A_B ( mxArray * A_m , mxArray * B_m )
2010-10-18 17:28:21 +02:00
{
unsigned int n_A = mxGetN ( A_m ) ;
unsigned int m_A = mxGetM ( A_m ) ;
2011-02-04 16:53:12 +01:00
double * A_d = mxGetPr ( A_m ) ;
2010-10-18 17:28:21 +02:00
unsigned int n_B = mxGetN ( B_m ) ;
2011-02-04 16:53:12 +01:00
double * B_d = mxGetPr ( B_m ) ;
mxArray * C_m = mxCreateDoubleMatrix ( m_A , n_B , mxREAL ) ;
double * C_d = mxGetPr ( C_m ) ;
2010-10-18 17:28:21 +02:00
for ( unsigned int j = 0 ; j < n_A ; j + + )
for ( unsigned int i = 0 ; i < m_A ; i + + )
{
unsigned int index = j * m_A + i ;
C_d [ index ] = A_d [ index ] - B_d [ index ] ;
2011-02-04 16:53:12 +01:00
}
2010-10-18 17:28:21 +02:00
return C_m ;
}
2011-02-04 16:53:12 +01:00
mxArray *
SparseMatrix : : Sparse_substract_A_SB ( mxArray * A_m , mxArray * B_m )
2010-10-18 17:28:21 +02:00
{
unsigned int n_B = mxGetN ( B_m ) ;
unsigned int m_B = mxGetM ( B_m ) ;
2011-02-04 16:53:12 +01:00
mwIndex * B_i = mxGetIr ( B_m ) ;
mwIndex * B_j = mxGetJc ( B_m ) ;
2010-10-18 17:28:21 +02:00
unsigned int total_nze_B = B_j [ n_B ] ;
2011-02-04 16:53:12 +01:00
double * B_d = mxGetPr ( B_m ) ;
mxArray * C_m = mxDuplicateArray ( A_m ) ;
double * C_d = mxGetPr ( C_m ) ;
2010-10-18 17:28:21 +02:00
unsigned int nze_B = 0 ;
unsigned int B_col = 0 ;
while ( nze_B < total_nze_B )
{
2011-02-04 16:53:12 +01:00
while ( nze_B > = ( unsigned int ) B_j [ B_col + 1 ] & & ( nze_B < total_nze_B ) )
2010-10-18 17:28:21 +02:00
B_col + + ;
C_d [ B_col * m_B + B_i [ nze_B ] ] - = B_d [ nze_B ] ;
nze_B + + ;
}
return C_m ;
}
2011-02-04 16:53:12 +01:00
mxArray *
SparseMatrix : : Sparse_substract_SA_SB ( mxArray * A_m , mxArray * B_m )
2010-10-18 17:28:21 +02:00
{
unsigned int n_A = mxGetN ( A_m ) ;
unsigned int m_A = mxGetM ( A_m ) ;
2011-02-04 16:53:12 +01:00
mwIndex * A_i = mxGetIr ( A_m ) ;
mwIndex * A_j = mxGetJc ( A_m ) ;
2010-10-18 17:28:21 +02:00
unsigned int total_nze_A = A_j [ n_A ] ;
2011-02-04 16:53:12 +01:00
double * A_d = mxGetPr ( A_m ) ;
2010-10-18 17:28:21 +02:00
unsigned int n_B = mxGetN ( B_m ) ;
2011-02-04 16:53:12 +01:00
mwIndex * B_i = mxGetIr ( B_m ) ;
mwIndex * B_j = mxGetJc ( B_m ) ;
2010-10-18 17:28:21 +02:00
unsigned int total_nze_B = B_j [ n_B ] ;
2011-02-04 16:53:12 +01:00
double * B_d = mxGetPr ( B_m ) ;
mxArray * C_m = mxCreateSparse ( m_A , n_B , m_A * n_B , mxREAL ) ;
mwIndex * C_i = mxGetIr ( C_m ) ;
mwIndex * C_j = mxGetJc ( C_m ) ;
double * C_d = mxGetPr ( C_m ) ;
2010-10-18 17:28:21 +02:00
unsigned int nze_B = 0 , nze_C = 0 , nze_A = 0 ;
unsigned int A_col = 0 , B_col = 0 , C_col = 0 ;
C_j [ C_col ] = 0 ;
while ( nze_A < total_nze_A | | nze_B < total_nze_B )
{
2011-02-04 16:53:12 +01:00
while ( nze_A > = ( unsigned int ) A_j [ A_col + 1 ] & & ( nze_A < total_nze_A ) )
2010-10-18 17:28:21 +02:00
A_col + + ;
int A_row = A_i [ nze_A ] ;
2011-02-04 16:53:12 +01:00
while ( nze_B > = ( unsigned int ) B_j [ B_col + 1 ] & & ( nze_B < total_nze_B ) )
2010-10-18 17:28:21 +02:00
B_col + + ;
int B_row = B_i [ nze_B ] ;
if ( A_col = = B_col )
{
if ( A_row = = B_row & & ( nze_B < total_nze_B & & nze_A < total_nze_A ) )
{
C_d [ nze_C ] = A_d [ nze_A + + ] - B_d [ nze_B + + ] ;
C_i [ nze_C ] = A_row ;
while ( C_col < A_col )
C_j [ + + C_col ] = nze_C ;
C_j [ A_col + 1 ] = nze_C + + ;
C_col = A_col ;
}
else if ( A_row < B_row | | ( nze_B > = total_nze_B & & nze_A < total_nze_A ) )
{
C_d [ nze_C ] = A_d [ nze_A + + ] ;
C_i [ nze_C ] = A_row ;
while ( C_col < A_col )
C_j [ + + C_col ] = nze_C ;
C_j [ A_col + 1 ] = nze_C + + ;
C_col = A_col ;
}
else
{
2011-02-04 16:53:12 +01:00
C_d [ nze_C ] = - B_d [ nze_B + + ] ;
2010-10-18 17:28:21 +02:00
C_i [ nze_C ] = B_row ;
while ( C_col < B_col )
C_j [ + + C_col ] = nze_C ;
C_j [ B_col + 1 ] = nze_C + + ;
C_col = B_col ;
}
}
else if ( A_col < B_col | | ( nze_B > = total_nze_B & & nze_A < total_nze_A ) )
{
C_d [ nze_C ] = A_d [ nze_A + + ] ;
C_i [ nze_C ] = A_row ;
while ( C_col < A_col )
C_j [ + + C_col ] = nze_C ;
C_j [ A_col + 1 ] = nze_C + + ;
C_col = A_col ;
}
else
{
2011-02-04 16:53:12 +01:00
C_d [ nze_C ] = - B_d [ nze_B + + ] ;
2010-10-18 17:28:21 +02:00
C_i [ nze_C ] = B_row ;
while ( C_col < B_col )
2011-02-04 16:53:12 +01:00
C_j [ + + C_col ] = nze_C ;
2010-10-18 17:28:21 +02:00
C_j [ B_col + 1 ] = nze_C + + ;
C_col = B_col ;
}
}
while ( C_col < n_B )
C_j [ + + C_col ] = nze_C ;
mxSetNzmax ( C_m , nze_C ) ;
return C_m ;
}
2011-02-04 16:53:12 +01:00
mxArray *
SparseMatrix : : mult_SAT_B ( mxArray * A_m , mxArray * B_m )
2010-10-18 17:28:21 +02:00
{
unsigned int n_A = mxGetN ( A_m ) ;
unsigned int m_A = mxGetM ( A_m ) ;
2011-02-04 16:53:12 +01:00
mwIndex * A_i = mxGetIr ( A_m ) ;
mwIndex * A_j = mxGetJc ( A_m ) ;
double * A_d = mxGetPr ( A_m ) ;
2010-10-18 17:28:21 +02:00
unsigned int n_B = mxGetN ( B_m ) ;
2011-02-04 16:53:12 +01:00
double * B_d = mxGetPr ( B_m ) ;
mxArray * C_m = mxCreateDoubleMatrix ( m_A , n_B , mxREAL ) ;
double * C_d = mxGetPr ( C_m ) ;
2010-10-18 17:28:21 +02:00
unsigned int nze_A = 0 ;
for ( unsigned int j = 0 ; j < n_B ; j + + )
{
for ( unsigned int i = 0 ; i < n_A ; i + + )
2011-02-04 16:53:12 +01:00
{
2010-10-18 17:28:21 +02:00
double sum = 0 ;
nze_A = A_j [ i ] ;
2011-02-04 16:53:12 +01:00
while ( nze_A < ( unsigned int ) A_j [ i + 1 ] )
2010-10-18 17:28:21 +02:00
{
unsigned int i_A = A_i [ nze_A ] ;
sum + = A_d [ nze_A + + ] * B_d [ i_A ] ;
}
C_d [ j * n_A + i ] = sum ;
}
}
return C_m ;
}
2011-02-04 16:53:12 +01:00
mxArray *
SparseMatrix : : Sparse_mult_SAT_B ( mxArray * A_m , mxArray * B_m )
2010-10-18 17:28:21 +02:00
{
unsigned int n_A = mxGetN ( A_m ) ;
unsigned int m_A = mxGetM ( A_m ) ;
2011-02-04 16:53:12 +01:00
mwIndex * A_i = mxGetIr ( A_m ) ;
mwIndex * A_j = mxGetJc ( A_m ) ;
double * A_d = mxGetPr ( A_m ) ;
2010-10-18 17:28:21 +02:00
unsigned int n_B = mxGetN ( B_m ) ;
unsigned int m_B = mxGetM ( B_m ) ;
2011-02-04 16:53:12 +01:00
double * B_d = mxGetPr ( B_m ) ;
mxArray * C_m = mxCreateSparse ( m_A , n_B , m_A * n_B , mxREAL ) ;
mwIndex * C_i = mxGetIr ( C_m ) ;
mwIndex * C_j = mxGetJc ( C_m ) ;
double * C_d = mxGetPr ( C_m ) ;
2010-10-18 17:28:21 +02:00
unsigned int nze_C = 0 , nze_A = 0 ;
unsigned int C_col = 0 ;
C_j [ C_col ] = 0 ;
for ( unsigned int j = 0 ; j < n_B ; j + + )
{
for ( unsigned int i = 0 ; i < n_A ; i + + )
2011-02-04 16:53:12 +01:00
{
2010-10-18 17:28:21 +02:00
double sum = 0 ;
nze_A = A_j [ i ] ;
2011-02-04 16:53:12 +01:00
while ( nze_A < ( unsigned int ) A_j [ i + 1 ] )
2010-10-18 17:28:21 +02:00
{
unsigned int i_A = A_i [ nze_A ] ;
sum + = A_d [ nze_A + + ] * B_d [ i_A ] ;
}
if ( fabs ( sum ) > 1e-10 )
{
C_d [ nze_C ] = sum ;
C_i [ nze_C ] = i ;
while ( C_col < j )
C_j [ + + C_col ] = nze_C ;
nze_C + + ;
}
}
}
while ( C_col < m_B )
C_j [ + + C_col ] = nze_C ;
mxSetNzmax ( C_m , nze_C ) ;
return C_m ;
}
2011-02-04 16:53:12 +01:00
mxArray *
SparseMatrix : : Sparse_mult_SAT_SB ( mxArray * A_m , mxArray * B_m )
2010-10-18 17:28:21 +02:00
{
unsigned int n_A = mxGetN ( A_m ) ;
unsigned int m_A = mxGetM ( A_m ) ;
2011-02-04 16:53:12 +01:00
mwIndex * A_i = mxGetIr ( A_m ) ;
mwIndex * A_j = mxGetJc ( A_m ) ;
double * A_d = mxGetPr ( A_m ) ;
2010-10-18 17:28:21 +02:00
unsigned int n_B = mxGetN ( B_m ) ;
2011-02-04 16:53:12 +01:00
mwIndex * B_i = mxGetIr ( B_m ) ;
mwIndex * B_j = mxGetJc ( B_m ) ;
double * B_d = mxGetPr ( B_m ) ;
mxArray * C_m = mxCreateSparse ( m_A , n_B , m_A * n_B , mxREAL ) ;
mwIndex * C_i = mxGetIr ( C_m ) ;
mwIndex * C_j = mxGetJc ( C_m ) ;
double * C_d = mxGetPr ( C_m ) ;
2010-10-18 17:28:21 +02:00
unsigned int nze_B = 0 , nze_C = 0 , nze_A = 0 ;
unsigned int C_col = 0 ;
C_j [ C_col ] = 0 ;
for ( unsigned int j = 0 ; j < n_B ; j + + )
{
for ( unsigned int i = 0 ; i < n_A ; i + + )
2011-02-04 16:53:12 +01:00
{
2010-10-18 17:28:21 +02:00
double sum = 0 ;
nze_B = B_j [ j ] ;
nze_A = A_j [ i ] ;
2011-02-04 16:53:12 +01:00
while ( nze_A < ( unsigned int ) A_j [ i + 1 ] & & nze_B < ( unsigned int ) B_j [ j + 1 ] )
2010-10-18 17:28:21 +02:00
{
unsigned int i_A = A_i [ nze_A ] ;
unsigned int i_B = B_i [ nze_B ] ;
2011-02-04 16:53:12 +01:00
if ( i_A = = i_B )
2010-10-18 17:28:21 +02:00
sum + = A_d [ nze_A + + ] * B_d [ nze_B + + ] ;
else if ( i_A < i_B )
nze_A + + ;
else
nze_B + + ;
}
if ( fabs ( sum ) > 1e-10 )
{
C_d [ nze_C ] = sum ;
C_i [ nze_C ] = i ;
while ( C_col < j )
C_j [ + + C_col ] = nze_C ;
nze_C + + ;
}
}
}
while ( C_col < n_B )
C_j [ + + C_col ] = nze_C ;
mxSetNzmax ( C_m , nze_C ) ;
return C_m ;
}
2011-02-04 16:53:12 +01:00
mxArray *
SparseMatrix : : Sparse_transpose ( mxArray * A_m )
2010-10-18 17:28:21 +02:00
{
unsigned int n_A = mxGetN ( A_m ) ;
unsigned int m_A = mxGetM ( A_m ) ;
2011-02-04 16:53:12 +01:00
mwIndex * A_i = mxGetIr ( A_m ) ;
mwIndex * A_j = mxGetJc ( A_m ) ;
2010-10-18 17:28:21 +02:00
unsigned int total_nze_A = A_j [ n_A ] ;
2011-02-04 16:53:12 +01:00
double * A_d = mxGetPr ( A_m ) ;
mxArray * C_m = mxCreateSparse ( n_A , m_A , total_nze_A , mxREAL ) ;
mwIndex * C_i = mxGetIr ( C_m ) ;
mwIndex * C_j = mxGetJc ( C_m ) ;
double * C_d = mxGetPr ( C_m ) ;
2010-10-18 17:28:21 +02:00
unsigned int nze_C = 0 , nze_A = 0 ;
memset ( C_j , 0 , m_A ) ;
map < pair < unsigned int , unsigned int > , double > B2 ;
for ( unsigned int i = 0 ; i < n_A ; i + + )
{
2011-02-04 16:53:12 +01:00
while ( nze_A < ( unsigned int ) A_j [ i + 1 ] )
2010-10-18 17:28:21 +02:00
{
C_j [ A_i [ nze_A ] + 1 ] + + ;
B2 [ make_pair ( A_i [ nze_A ] , i ) ] = A_d [ nze_A ] ;
nze_A + + ;
}
}
for ( unsigned int i = 0 ; i < m_A ; i + + )
C_j [ i + 1 ] + = C_j [ i ] ;
for ( map < pair < unsigned int , unsigned int > , double > : : const_iterator it = B2 . begin ( ) ; it ! = B2 . end ( ) ; it + + )
{
C_d [ nze_C ] = it - > second ;
C_i [ nze_C + + ] = it - > first . second ;
}
return C_m ;
}
2010-09-24 12:52:58 +02:00
void
2011-02-04 16:53:12 +01:00
SparseMatrix : : Solve_Matlab_Relaxation ( mxArray * A_m , mxArray * b_m , unsigned int Size , double slowc_l , bool is_two_boundaries , int it_ )
2010-09-24 12:52:58 +02:00
{
mxArray * B1 , * C1 , * A2 , * B2 , * A3 , * b1 , * b2 ;
2011-02-04 16:53:12 +01:00
double * b_m_d = mxGetPr ( b_m ) ;
2010-09-24 12:52:58 +02:00
if ( ! b_m_d )
2009-01-20 23:04:37 +01:00
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Solve_Matlab_Relaxation, can't retrieve b_m vector \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2010-09-17 09:57:38 +02:00
}
mwIndex * A_m_i = mxGetIr ( A_m ) ;
if ( ! A_m_i )
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Solve_Matlab_Relaxation, can't allocate A_m_i index vector \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2010-09-17 09:57:38 +02:00
}
mwIndex * A_m_j = mxGetJc ( A_m ) ;
if ( ! A_m_j )
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Solve_Matlab_Relaxation, can't allocate A_m_j index vector \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2010-09-17 09:57:38 +02:00
}
double * A_m_d = mxGetPr ( A_m ) ;
if ( ! A_m_d )
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Solve_Matlab_Relaxation, can't retrieve A matrix \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2010-09-17 09:57:38 +02:00
}
unsigned int max_nze = A_m_j [ Size * periods ] ;
unsigned int nze = 0 ;
unsigned int var = A_m_j [ nze ] ;
B1 = mxCreateSparse ( Size , Size , Size * Size , mxREAL ) ;
2011-02-04 16:53:12 +01:00
mwIndex * B1_i = mxGetIr ( B1 ) ;
mwIndex * B1_j = mxGetJc ( B1 ) ;
double * B1_d = mxGetPr ( B1 ) ;
2010-09-17 09:57:38 +02:00
unsigned int B1_nze = 0 ;
unsigned int B1_var = 0 ;
B1_i [ B1_nze ] = 0 ;
B1_j [ B1_var ] = 0 ;
C1 = mxCreateSparse ( Size , Size , Size * Size , mxREAL ) ;
2011-02-04 16:53:12 +01:00
mwIndex * C1_i = mxGetIr ( C1 ) ;
mwIndex * C1_j = mxGetJc ( C1 ) ;
double * C1_d = mxGetPr ( C1 ) ;
2010-09-17 09:57:38 +02:00
unsigned int C1_nze = 0 ;
unsigned int C1_var = 0 ;
C1_i [ C1_nze ] = 0 ;
C1_j [ C1_var ] = 0 ;
A2 = mxCreateSparse ( Size , Size , Size * Size , mxREAL ) ;
2011-02-04 16:53:12 +01:00
mwIndex * A2_i = mxGetIr ( A2 ) ;
mwIndex * A2_j = mxGetJc ( A2 ) ;
double * A2_d = mxGetPr ( A2 ) ;
2010-09-17 09:57:38 +02:00
unsigned int A2_nze = 0 ;
unsigned int A2_var = 0 ;
A2_i [ A2_nze ] = 0 ;
A2_j [ A2_var ] = 0 ;
B2 = mxCreateSparse ( Size , Size , Size * Size , mxREAL ) ;
2011-02-04 16:53:12 +01:00
mwIndex * B2_i = mxGetIr ( B2 ) ;
mwIndex * B2_j = mxGetJc ( B2 ) ;
double * B2_d = mxGetPr ( B2 ) ;
2010-09-17 09:57:38 +02:00
unsigned int B2_nze = 0 ;
unsigned int B2_var = 0 ;
B2_i [ B2_nze ] = 0 ;
B2_j [ B2_var ] = 0 ;
A3 = mxCreateSparse ( Size , Size , Size * Size , mxREAL ) ;
2011-02-04 16:53:12 +01:00
mwIndex * A3_i = mxGetIr ( A3 ) ;
mwIndex * A3_j = mxGetJc ( A3 ) ;
double * A3_d = mxGetPr ( A3 ) ;
2010-09-17 09:57:38 +02:00
unsigned int A3_nze = 0 ;
unsigned int A3_var = 0 ;
A3_i [ A3_nze ] = 0 ;
A3_j [ A3_var ] = 0 ;
b1 = mxCreateDoubleMatrix ( Size , 1 , mxREAL ) ;
2011-02-04 16:53:12 +01:00
double * b1_d = mxGetPr ( b1 ) ;
2010-09-17 09:57:38 +02:00
b2 = mxCreateDoubleMatrix ( Size , 1 , mxREAL ) ;
2011-02-04 16:53:12 +01:00
double * b2_d = mxGetPr ( b2 ) ;
2010-09-17 09:57:38 +02:00
unsigned int eq = 0 ;
/*B1 C1
A2 B2
2011-02-04 16:53:12 +01:00
A3 */
2010-09-17 09:57:38 +02:00
while ( var < 2 * Size & & nze < max_nze )
{
2011-02-04 16:53:12 +01:00
if ( ( unsigned int ) A_m_j [ var + 1 ] < = nze )
2010-09-17 09:57:38 +02:00
{
if ( var < Size )
b1_d [ var ] = b_m_d [ var ] ;
else
b2_d [ var - Size ] = b_m_d [ var ] ;
var + + ;
}
eq = A_m_i [ nze ] ;
if ( var < Size )
{
if ( eq < Size )
{
while ( B1_var < var )
B1_j [ + + B1_var ] = B1_nze ;
B1_i [ B1_nze ] = eq ;
B1_d [ B1_nze ] = A_m_d [ nze ] ;
B1_nze + + ;
}
else
{
while ( A2_var < var )
A2_j [ + + A2_var ] = A2_nze ;
A2_i [ A2_nze ] = eq - Size ;
A2_d [ A2_nze ] = A_m_d [ nze ] ;
A2_nze + + ;
}
}
2011-02-04 16:53:12 +01:00
else if ( var < 2 * Size )
2010-09-17 09:57:38 +02:00
{
if ( eq < Size )
{
while ( C1_var < var - Size )
C1_j [ + + C1_var ] = C1_nze ;
C1_i [ C1_nze ] = eq ;
C1_d [ C1_nze ] = A_m_d [ nze ] ;
C1_nze + + ;
}
else if ( eq < 2 * Size )
{
while ( B2_var < var - Size )
B2_j [ + + B2_var ] = B2_nze ;
B2_i [ B2_nze ] = eq - Size ;
B2_d [ B2_nze ] = A_m_d [ nze ] ;
B2_nze + + ;
}
else
{
while ( A3_var < var - Size )
A3_j [ + + A3_var ] = A3_nze ;
A3_i [ A3_nze ] = eq - 2 * Size ;
A3_d [ A3_nze ] = A_m_d [ nze ] ;
A3_nze + + ;
}
}
nze + + ;
}
while ( B1_var < Size )
B1_j [ + + B1_var ] = B1_nze ;
while ( C1_var < Size )
C1_j [ + + C1_var ] = C1_nze ;
while ( A2_var < Size )
A2_j [ + + A2_var ] = A2_nze ;
while ( B2_var < Size )
B2_j [ + + B2_var ] = B2_nze ;
2011-02-04 16:53:12 +01:00
while ( A3_var < Size )
2010-09-17 09:57:38 +02:00
A3_j [ + + A3_var ] = A3_nze ;
2011-01-14 19:22:29 +01:00
mxArray * d1 = NULL ;
2011-02-04 16:53:12 +01:00
vector < pair < mxArray * , mxArray * > > triangular_form ;
double sumc = 0 , C_sumc = 1000 ;
2010-10-18 17:28:21 +02:00
mxArray * B1_inv = NULL ;
mxArray * B1_inv_t = NULL ;
2010-09-17 09:57:38 +02:00
for ( int t = 1 ; t < = periods ; t + + )
{
if ( abs ( sumc / C_sumc - 1 ) > 1e-10 * res1 )
{
C_sumc = sumc ;
2010-10-18 17:28:21 +02:00
if ( B1_inv )
mxDestroyArray ( B1_inv ) ;
2010-09-17 09:57:38 +02:00
mexCallMATLAB ( 1 , & B1_inv , 1 , & B1 , " inv " ) ;
mwIndex * B_inv_j = mxGetJc ( B1_inv ) ;
unsigned int B_inv_nze = B_inv_j [ Size ] ;
2011-02-04 16:53:12 +01:00
double * B_inv_d = mxGetPr ( B1_inv ) ;
2010-09-17 09:57:38 +02:00
sumc = 0 ;
2011-02-04 16:53:12 +01:00
for ( unsigned int i = 0 ; i < B_inv_nze ; i + + )
2010-09-17 09:57:38 +02:00
sumc + = fabs ( B_inv_d [ i ] ) ;
}
2010-10-18 17:28:21 +02:00
B1_inv_t = Sparse_transpose ( B1_inv ) ;
mxArray * S1 = Sparse_mult_SAT_SB ( B1_inv_t , C1 ) ;
2010-09-17 09:57:38 +02:00
2010-10-18 17:28:21 +02:00
d1 = mult_SAT_B ( B1_inv_t , b1 ) ;
2010-09-17 09:57:38 +02:00
if ( t < periods )
2011-02-04 16:53:12 +01:00
//Computation for the next lines
2010-09-17 09:57:38 +02:00
{
2010-10-18 17:28:21 +02:00
mxDestroyArray ( B1_inv_t ) ;
mxArray * A2_t = Sparse_transpose ( A2 ) ;
mxDestroyArray ( A2 ) ;
mxArray * tmp = Sparse_mult_SAT_SB ( A2_t , S1 ) ;
2010-09-17 09:57:38 +02:00
mxDestroyArray ( B1 ) ;
2010-10-18 17:28:21 +02:00
B1 = Sparse_substract_SA_SB ( B2 , tmp ) ;
2010-09-17 09:57:38 +02:00
mxDestroyArray ( tmp ) ;
2010-10-18 17:28:21 +02:00
tmp = mult_SAT_B ( A2_t , d1 ) ;
b1 = substract_A_B ( b2 , tmp ) ;
2010-09-17 09:57:38 +02:00
mxDestroyArray ( tmp ) ;
2010-10-18 17:28:21 +02:00
triangular_form . push_back ( make_pair ( S1 , d1 ) ) ;
mxDestroyArray ( A2_t ) ;
2010-09-17 09:57:38 +02:00
}
A2 = mxDuplicateArray ( A3 ) ;
2010-10-18 17:28:21 +02:00
2010-09-17 09:57:38 +02:00
//I S1
//0 B1 C1 =>B1 =
// A2 B2 => A2 = A3
// A3
C1_nze = B2_nze = A3_nze = 0 ;
C1_var = B2_var = A3_var = 0 ;
if ( nze < max_nze )
nze - - ;
while ( var < ( t + 2 ) * Size & & nze < max_nze )
{
2011-02-04 16:53:12 +01:00
if ( ( unsigned int ) A_m_j [ var + 1 ] < = nze )
2010-09-17 09:57:38 +02:00
{
b2_d [ var - ( t + 1 ) * Size ] = b_m_d [ var ] ;
var + + ;
}
eq = A_m_i [ nze ] ;
if ( eq < ( t + 1 ) * Size )
{
C1_d [ C1_nze ] = A_m_d [ nze ] ;
C1_nze + + ;
}
else if ( eq < ( t + 2 ) * Size )
{
B2_d [ B2_nze ] = A_m_d [ nze ] ;
B2_nze + + ;
}
else
{
A3_d [ A3_nze ] = A_m_d [ nze ] ;
A3_nze + + ;
}
nze + + ;
}
}
double * d1_d = mxGetPr ( d1 ) ;
2011-02-04 16:53:12 +01:00
for ( unsigned i = 0 ; i < Size ; i + + )
2010-09-17 09:57:38 +02:00
{
int eq = index_vara [ i + Size * ( y_kmin + periods - 1 ) ] ;
2011-02-04 16:53:12 +01:00
double yy = - ( d1_d [ i ] + y [ eq ] ) ;
2010-09-17 09:57:38 +02:00
direction [ eq ] = yy ;
y [ eq ] + = slowc_l * yy ;
}
2011-02-04 16:53:12 +01:00
pair < mxArray * , mxArray * > tf ;
2010-09-17 09:57:38 +02:00
for ( int t = periods - 2 ; t > = 0 ; t - - )
{
2011-02-04 16:53:12 +01:00
mxArray * tmp ;
2010-09-17 09:57:38 +02:00
tf = triangular_form . back ( ) ;
triangular_form . pop_back ( ) ;
2011-02-04 16:53:12 +01:00
mxArray * tf_first_t = Sparse_transpose ( tf . first ) ;
2010-10-18 17:28:21 +02:00
mxDestroyArray ( tf . first ) ;
tmp = mult_SAT_B ( tf_first_t , d1 ) ;
d1 = substract_A_B ( tf . second , tmp ) ;
2010-09-17 09:57:38 +02:00
d1_d = mxGetPr ( d1 ) ;
mxDestroyArray ( tmp ) ;
2011-02-04 16:53:12 +01:00
for ( unsigned i = 0 ; i < Size ; i + + )
2010-09-17 09:57:38 +02:00
{
int eq = index_vara [ i + Size * ( y_kmin + t ) ] ;
2011-02-04 16:53:12 +01:00
double yy = - ( d1_d [ i ] + y [ eq ] ) ;
2010-09-17 09:57:38 +02:00
direction [ eq ] = yy ;
y [ eq ] + = slowc_l * yy ;
}
2010-10-18 17:28:21 +02:00
mxDestroyArray ( tf_first_t ) ;
2010-09-17 09:57:38 +02:00
mxDestroyArray ( tf . second ) ;
}
mxDestroyArray ( B1 ) ;
mxDestroyArray ( C1 ) ;
mxDestroyArray ( A2 ) ;
mxDestroyArray ( B2 ) ;
mxDestroyArray ( A3 ) ;
mxDestroyArray ( b1 ) ;
mxDestroyArray ( b2 ) ;
mxDestroyArray ( A_m ) ;
mxDestroyArray ( b_m ) ;
}
2010-07-23 11:20:24 +02:00
void
2011-02-04 16:53:12 +01:00
SparseMatrix : : Solve_Matlab_LU_UMFPack ( mxArray * A_m , mxArray * b_m , int Size , double slowc_l , bool is_two_boundaries , int it_ )
2010-07-23 11:20:24 +02:00
{
int n = mxGetM ( A_m ) ;
mxArray * z ;
mxArray * rhs [ 2 ] ;
rhs [ 0 ] = A_m ;
rhs [ 1 ] = b_m ;
2010-08-18 13:51:57 +02:00
mexCallMATLAB ( 1 , & z , 2 , rhs , " mldivide " ) ;
2010-07-23 11:20:24 +02:00
double * res = mxGetPr ( z ) ;
2010-08-18 13:51:57 +02:00
if ( is_two_boundaries )
for ( int i = 0 ; i < n ; i + + )
{
int eq = index_vara [ i + Size * y_kmin ] ;
2011-02-04 16:53:12 +01:00
double yy = - ( res [ i ] + y [ eq ] ) ;
direction [ eq ] = yy ;
2010-08-18 13:51:57 +02:00
y [ eq ] + = slowc_l * yy ;
}
else
for ( int i = 0 ; i < n ; i + + )
{
2010-09-17 09:57:38 +02:00
int eq = index_vara [ i ] ;
2011-02-04 16:53:12 +01:00
double yy = - ( res [ i ] + y [ eq + it_ * y_size ] ) ;
2010-08-18 13:51:57 +02:00
direction [ eq ] = yy ;
y [ eq + it_ * y_size ] + = slowc_l * yy ;
}
2010-07-23 11:20:24 +02:00
mxDestroyArray ( A_m ) ;
mxDestroyArray ( b_m ) ;
mxDestroyArray ( z ) ;
}
void
2011-02-04 16:53:12 +01:00
SparseMatrix : : Solve_Matlab_GMRES ( mxArray * A_m , mxArray * b_m , int Size , double slowc , int block , bool is_two_boundaries , int it_ , bool steady_state , mxArray * x0_m )
2010-07-23 11:20:24 +02:00
{
2010-10-18 17:28:21 +02:00
# ifdef OCTAVE_MEX_FILE
ostringstream tmp ;
if ( steady_state )
2010-10-22 16:49:47 +02:00
tmp < < " GMRES method is not implemented in Octave. You cannot use solve_algo=7, change solve_algo. \n " ;
2010-10-18 17:28:21 +02:00
else
tmp < < " GMRES method is not implemented in Octave. You cannot use stack_solve_algo=2, change stack_solve_algo. \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
# endif
2010-07-23 11:20:24 +02:00
int n = mxGetM ( A_m ) ;
mxArray * lhs0 [ 2 ] ;
mxArray * rhs0 [ 2 ] ;
rhs0 [ 0 ] = A_m ;
rhs0 [ 1 ] = mxCreateDoubleScalar ( lu_inc_tol ) ;
mexCallMATLAB ( 2 , lhs0 , 2 , rhs0 , " luinc " ) ;
mxArray * L1 = lhs0 [ 0 ] ;
mxArray * U1 = lhs0 [ 1 ] ;
/*[za,flag1] = gmres(g1a,b,Blck_size,1e-6,Blck_size*periods,L1,U1);*/
2011-01-03 14:15:12 +01:00
mxArray * rhs [ 8 ] ;
2010-07-23 11:20:24 +02:00
rhs [ 0 ] = A_m ;
rhs [ 1 ] = b_m ;
rhs [ 2 ] = mxCreateDoubleScalar ( Size ) ;
rhs [ 3 ] = mxCreateDoubleScalar ( 1e-6 ) ;
rhs [ 4 ] = mxCreateDoubleScalar ( n ) ;
rhs [ 5 ] = L1 ;
rhs [ 6 ] = U1 ;
2011-01-03 14:15:12 +01:00
rhs [ 7 ] = x0_m ;
2010-07-23 11:20:24 +02:00
mxArray * lhs [ 2 ] ;
2011-02-04 16:53:12 +01:00
mexCallMATLAB ( 2 , lhs , 8 , rhs , " gmres " ) ;
2010-07-23 11:20:24 +02:00
mxArray * z = lhs [ 0 ] ;
mxArray * flag = lhs [ 1 ] ;
double * flag1 = mxGetPr ( flag ) ;
mxDestroyArray ( rhs0 [ 1 ] ) ;
mxDestroyArray ( rhs [ 2 ] ) ;
mxDestroyArray ( rhs [ 3 ] ) ;
mxDestroyArray ( rhs [ 4 ] ) ;
mxDestroyArray ( rhs [ 5 ] ) ;
mxDestroyArray ( rhs [ 6 ] ) ;
2010-10-18 17:28:21 +02:00
if ( * flag1 > 0 )
2010-07-23 11:20:24 +02:00
{
ostringstream tmp ;
if ( * flag1 = = 1 )
{
2010-10-18 17:28:21 +02:00
tmp < < " Error in bytecode: No convergence inside GMRES, in block " < < block + 1 ;
2010-07-23 11:20:24 +02:00
mexWarnMsgTxt ( tmp . str ( ) . c_str ( ) ) ;
}
else if ( * flag1 = = 2 )
{
2010-10-18 17:28:21 +02:00
tmp < < " Error in bytecode: Preconditioner is ill-conditioned, in block " < < block + 1 ;
2010-07-23 11:20:24 +02:00
mexWarnMsgTxt ( tmp . str ( ) . c_str ( ) ) ;
}
else if ( * flag1 = = 3 )
{
2010-10-18 17:28:21 +02:00
tmp < < " Error in bytecode: GMRES stagnated (Two consecutive iterates were the same.), in block " < < block + 1 ;
2010-07-23 11:20:24 +02:00
mexWarnMsgTxt ( tmp . str ( ) . c_str ( ) ) ;
}
lu_inc_tol / = 10 ;
}
else
{
double * res = mxGetPr ( z ) ;
2010-08-18 13:51:57 +02:00
if ( is_two_boundaries )
for ( int i = 0 ; i < n ; i + + )
{
int eq = index_vara [ i + Size * y_kmin ] ;
2011-02-04 16:53:12 +01:00
double yy = - ( res [ i ] + y [ eq ] ) ;
2010-08-18 13:51:57 +02:00
direction [ eq ] = yy ;
y [ eq ] + = slowc * yy ;
}
else
for ( int i = 0 ; i < n ; i + + )
{
int eq = index_vara [ i ] ;
2011-02-04 16:53:12 +01:00
double yy = - ( res [ i ] + y [ eq + it_ * y_size ] ) ;
2010-10-18 17:28:21 +02:00
direction [ eq ] = yy ;
y [ eq + it_ * y_size ] + = slowc * yy ;
2010-08-18 13:51:57 +02:00
}
2010-07-23 11:20:24 +02:00
}
mxDestroyArray ( A_m ) ;
mxDestroyArray ( b_m ) ;
mxDestroyArray ( z ) ;
mxDestroyArray ( flag ) ;
}
void
2011-02-04 16:53:12 +01:00
SparseMatrix : : Solve_Matlab_BiCGStab ( mxArray * A_m , mxArray * b_m , int Size , double slowc , int block , bool is_two_boundaries , int it_ , mxArray * x0_m , bool steady_state )
2010-07-23 11:20:24 +02:00
{
2011-01-14 19:22:29 +01:00
unsigned int n = mxGetM ( A_m ) ;
2010-07-23 11:20:24 +02:00
/*[L1, U1]=luinc(g1a,luinc_tol);*/
mxArray * lhs0 [ 2 ] ;
mxArray * rhs0 [ 2 ] ;
rhs0 [ 0 ] = A_m ;
rhs0 [ 1 ] = mxCreateDoubleScalar ( lu_inc_tol ) ;
mexCallMATLAB ( 2 , lhs0 , 2 , rhs0 , " luinc " ) ;
mxArray * L1 = lhs0 [ 0 ] ;
mxArray * U1 = lhs0 [ 1 ] ;
2011-01-03 14:15:12 +01:00
double flags = 2 ;
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*/
{
2011-02-04 16:53:12 +01:00
mxArray * res = mult_SAT_B ( Sparse_transpose ( A_m ) , x0_m ) ;
2011-01-03 14:15:12 +01:00
double * resid = mxGetPr ( res ) ;
double * b = mxGetPr ( b_m ) ;
for ( unsigned int i = 0 ; i < n ; i + + )
resid [ i ] = b [ i ] - resid [ i ] ;
mxArray * rhs [ 2 ] ;
mxArray * lhs [ 1 ] ;
rhs [ 0 ] = L1 ;
rhs [ 1 ] = res ;
mexCallMATLAB ( 1 , lhs , 2 , rhs , " mldivide " ) ;
rhs [ 0 ] = U1 ;
rhs [ 1 ] = lhs [ 0 ] ;
mexCallMATLAB ( 1 , lhs , 2 , rhs , " mldivide " ) ;
z = lhs [ 0 ] ;
double * phat = mxGetPr ( z ) ;
double * x0 = mxGetPr ( x0_m ) ;
for ( unsigned int i = 0 ; i < n ; i + + )
phat [ i ] = x0 [ i ] + phat [ i ] ;
/*Check the solution*/
res = mult_SAT_B ( Sparse_transpose ( A_m ) , z ) ;
resid = mxGetPr ( res ) ;
double cum_abs = 0 ;
for ( unsigned int i = 0 ; i < n ; i + + )
{
resid [ i ] = b [ i ] - resid [ i ] ;
cum_abs + = fabs ( resid [ i ] ) ;
}
//mexPrintf("cum_abs=%g\n", cum_abs);
if ( cum_abs > 1e-7 )
flags = 2 ;
else
flags = 0 ;
mxDestroyArray ( res ) ;
}
//else
if ( flags = = 2 )
{
/*[za,flag1] = bicgstab(g1a,b,1e-6,Blck_size*periods,L1,U1);*/
mxArray * rhs [ 7 ] ;
rhs [ 0 ] = A_m ;
rhs [ 1 ] = b_m ;
rhs [ 2 ] = mxCreateDoubleScalar ( 1e-6 ) ;
rhs [ 3 ] = mxCreateDoubleScalar ( n ) ;
rhs [ 4 ] = L1 ;
rhs [ 5 ] = U1 ;
rhs [ 6 ] = x0_m ;
mxArray * lhs [ 2 ] ;
2011-02-04 16:53:12 +01:00
mexCallMATLAB ( 2 , lhs , 7 , rhs , " bicgstab " ) ;
2011-01-03 14:15:12 +01:00
z = lhs [ 0 ] ;
mxArray * flag = lhs [ 1 ] ;
double * flag1 = mxGetPr ( flag ) ;
flags = flag1 [ 0 ] ;
mxDestroyArray ( flag ) ;
mxDestroyArray ( rhs [ 2 ] ) ;
mxDestroyArray ( rhs [ 3 ] ) ;
mxDestroyArray ( rhs [ 4 ] ) ;
mxDestroyArray ( rhs [ 5 ] ) ;
}
/*mexPrintf("z");
2011-02-04 16:53:12 +01:00
mexCallMATLAB ( 0 , NULL , 1 , & z , " disp " ) ; */
2010-07-23 11:20:24 +02:00
mxDestroyArray ( rhs0 [ 1 ] ) ;
2011-01-03 14:15:12 +01:00
if ( flags > 0 )
2010-07-23 11:20:24 +02:00
{
ostringstream tmp ;
2011-01-03 14:15:12 +01:00
if ( flags = = 1 )
2010-07-23 11:20:24 +02:00
{
2010-10-18 17:28:21 +02:00
tmp < < " Error in bytecode: No convergence inside BiCGStab, in block " < < block + 1 ;
2010-07-23 11:20:24 +02:00
mexWarnMsgTxt ( tmp . str ( ) . c_str ( ) ) ;
}
2011-01-03 14:15:12 +01:00
else if ( flags = = 2 )
2010-07-23 11:20:24 +02:00
{
2010-10-18 17:28:21 +02:00
tmp < < " Error in bytecode: Preconditioner is ill-conditioned, in block " < < block + 1 ;
2010-07-23 11:20:24 +02:00
mexWarnMsgTxt ( tmp . str ( ) . c_str ( ) ) ;
}
2011-01-03 14:15:12 +01:00
else if ( flags = = 3 )
2010-07-23 11:20:24 +02:00
{
2010-10-18 17:28:21 +02:00
tmp < < " Error in bytecode: BiCGStab stagnated (Two consecutive iterates were the same.), in block " < < block + 1 ;
2010-07-23 11:20:24 +02:00
mexWarnMsgTxt ( tmp . str ( ) . c_str ( ) ) ;
}
lu_inc_tol / = 10 ;
}
else
{
double * res = mxGetPr ( z ) ;
2010-08-18 13:51:57 +02:00
if ( is_two_boundaries )
2011-01-14 19:22:29 +01:00
for ( unsigned int i = 0 ; i < n ; i + + )
2010-08-18 13:51:57 +02:00
{
int eq = index_vara [ i + Size * y_kmin ] ;
2011-02-04 16:53:12 +01:00
double yy = - ( res [ i ] + y [ eq ] ) ;
direction [ eq ] = yy ;
2010-08-18 13:51:57 +02:00
y [ eq ] + = slowc * yy ;
}
else
2011-01-14 19:22:29 +01:00
for ( unsigned int i = 0 ; i < n ; i + + )
2010-08-18 13:51:57 +02:00
{
int eq = index_vara [ i ] ;
2011-02-04 16:53:12 +01:00
double yy = - ( res [ i ] + y [ eq + it_ * y_size ] ) ;
direction [ eq ] = yy ;
2010-10-18 17:28:21 +02:00
y [ eq + it_ * y_size ] + = slowc * yy ;
2010-08-18 13:51:57 +02:00
}
2010-07-23 11:20:24 +02:00
}
mxDestroyArray ( A_m ) ;
mxDestroyArray ( b_m ) ;
mxDestroyArray ( z ) ;
}
2010-09-24 12:52:58 +02:00
void
2012-03-06 11:03:13 +01:00
SparseMatrix : : Singular_display ( int block , int Size , bool steady_state , it_code_type it_code )
{
bool zero_solution ;
Simple_Init ( Size , IM_i , zero_solution ) ;
NonZeroElem * first ;
mxArray * rhs [ 1 ] ;
rhs [ 0 ] = mxCreateDoubleMatrix ( Size , Size , mxREAL ) ;
double * pind ;
pind = mxGetPr ( rhs [ 0 ] ) ;
for ( int j = 0 ; j < Size * Size ; j + + )
pind [ j ] = 0.0 ;
for ( int ii = 0 ; ii < Size ; ii + + )
{
int nb_eq = At_Col ( ii , & first ) ;
for ( int j = 0 ; j < nb_eq ; j + + )
{
int k = first - > u_index ;
int jj = first - > r_index ;
pind [ ii * Size + jj ] = u [ k ] ;
first = first - > NZE_C_N ;
}
}
mxArray * lhs [ 3 ] ;
mexCallMATLAB ( 3 , lhs , 1 , rhs , " svd " ) ;
mxArray * SVD_u = lhs [ 0 ] ;
mxArray * SVD_s = lhs [ 1 ] ;
2012-06-06 16:30:28 +02:00
//mxArray* SVD_v = lhs[2];
2012-03-06 11:03:13 +01:00
double * SVD_ps = mxGetPr ( SVD_s ) ;
double * SVD_pu = mxGetPr ( SVD_u ) ;
for ( int i = 0 ; i < Size ; i + + )
{
if ( abs ( SVD_ps [ i * ( 1 + Size ) ] ) < 1e-12 )
{
mexPrintf ( " The following equations form a linear combination: \n " ) ;
double max_u = 0 ;
for ( int j = 0 ; j < Size ; j + + )
if ( abs ( SVD_pu [ j + i * Size ] ) > abs ( max_u ) )
max_u = SVD_pu [ j + i * Size ] ;
vector < int > equ_list ;
for ( int j = 0 ; j < Size ; j + + )
{
double rr = SVD_pu [ j + i * Size ] / max_u ;
if ( rr < - 1e-10 )
{
equ_list . push_back ( j ) ;
if ( rr ! = - 1 )
mexPrintf ( " - %3.2f*Dequ_%d_dy " , abs ( rr ) , j + 1 ) ;
else
mexPrintf ( " - Dequ_%d_dy " , j + 1 ) ;
}
else if ( rr > 1e-10 )
{
equ_list . push_back ( j ) ;
if ( j > 0 )
if ( rr ! = 1 )
mexPrintf ( " + %3.2f*Dequ_%d_dy " , rr , j + 1 ) ;
else
mexPrintf ( " + Dequ_%d_dy " , j + 1 ) ;
else
if ( rr ! = 1 )
mexPrintf ( " %3.2f*Dequ_%d_dy " , rr , j + 1 ) ;
else
mexPrintf ( " Dequ_%d_dy " , j + 1 ) ;
}
}
mexPrintf ( " = 0 \n " ) ;
/*mexPrintf(" with:\n");
it_code = get_begin_block ( block ) ;
for ( int j = 0 ; j < Size ; j + + )
{
if ( find ( equ_list . begin ( ) , equ_list . end ( ) , j ) ! = equ_list . end ( ) )
mexPrintf ( " equ_%d: %s \n " , j , print_expression ( it_code_expr , false , Size , block , steady_state , 0 , 0 , it_code , true ) . c_str ( ) ) ;
} */
}
}
mxDestroyArray ( lhs [ 0 ] ) ;
mxDestroyArray ( lhs [ 1 ] ) ;
mxDestroyArray ( lhs [ 2 ] ) ;
ostringstream tmp ;
if ( block > 1 )
tmp < < " in Solve_ByteCode_Sparse_GaussianElimination, singular system in block " < < block + 1 < < " \n " ;
else
tmp < < " in Solve_ByteCode_Sparse_GaussianElimination, singular system \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
}
bool
2010-09-24 12:52:58 +02:00
SparseMatrix : : Solve_ByteCode_Sparse_GaussianElimination ( int Size , int blck , bool steady_state , int it_ )
2007-11-21 00:27:30 +01:00
{
2010-09-24 12:52:58 +02:00
bool one ;
2009-08-25 11:43:01 +02:00
int pivj = 0 , pivk = 0 ;
2010-09-24 12:52:58 +02:00
double * piv_v ;
int * pivj_v , * pivk_v , * NR ;
int l , N_max ;
NonZeroElem * first , * firsta , * first_suba ;
2009-08-25 11:43:01 +02:00
double piv_abs ;
2010-09-24 12:52:58 +02:00
NonZeroElem * * bc ;
bc = ( NonZeroElem * * ) mxMalloc ( Size * sizeof ( * bc ) ) ;
piv_v = ( double * ) mxMalloc ( Size * sizeof ( double ) ) ;
pivj_v = ( int * ) mxMalloc ( Size * sizeof ( int ) ) ;
pivk_v = ( int * ) mxMalloc ( Size * sizeof ( int ) ) ;
NR = ( int * ) mxMalloc ( Size * sizeof ( int ) ) ;
for ( int i = 0 ; i < Size ; i + + )
{
/*finding the max-pivot*/
double piv = piv_abs = 0 ;
int nb_eq = At_Col ( i , & first ) ;
l = 0 ; N_max = 0 ;
one = false ;
piv_abs = 0 ;
for ( int j = 0 ; j < nb_eq ; j + + )
{
if ( ! line_done [ first - > r_index ] )
{
int k = first - > u_index ;
int jj = first - > r_index ;
int NRow_jj = NRow ( jj ) ;
piv_v [ l ] = u [ k ] ;
double piv_fabs = fabs ( u [ k ] ) ;
pivj_v [ l ] = jj ;
pivk_v [ l ] = k ;
NR [ l ] = NRow_jj ;
if ( NRow_jj = = 1 & & ! one )
{
one = true ;
piv_abs = piv_fabs ;
N_max = NRow_jj ;
}
if ( ! one )
{
if ( piv_fabs > piv_abs )
piv_abs = piv_fabs ;
if ( NRow_jj > N_max )
N_max = NRow_jj ;
}
else
{
if ( NRow_jj = = 1 )
{
if ( piv_fabs > piv_abs )
piv_abs = piv_fabs ;
if ( NRow_jj > N_max )
N_max = NRow_jj ;
}
}
l + + ;
}
first = first - > NZE_C_N ;
}
if ( piv_abs < eps )
{
mxFree ( piv_v ) ;
mxFree ( pivj_v ) ;
mxFree ( pivk_v ) ;
mxFree ( NR ) ;
mxFree ( bc ) ;
if ( steady_state )
{
if ( blck > 1 )
mexPrintf ( " Error: singular system in Simulate_NG in block %d \n " , blck + 1 ) ;
else
mexPrintf ( " Error: singular system in Simulate_NG \n " ) ;
2012-03-06 11:03:13 +01:00
return true ;
2010-09-24 12:52:58 +02:00
}
else
{
ostringstream tmp ;
if ( blck > 1 )
tmp < < " in Solve_ByteCode_Sparse_GaussianElimination, singular system in block " < < blck + 1 < < " \n " ;
else
tmp < < " in Solve_ByteCode_Sparse_GaussianElimination, singular system \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
}
}
double markovitz = 0 , markovitz_max = - 9e70 ;
if ( ! one )
{
for ( int j = 0 ; j < l ; j + + )
{
if ( N_max > 0 & & NR [ j ] > 0 )
{
if ( fabs ( piv_v [ j ] ) > 0 )
{
if ( markowitz_c > 0 )
markovitz = exp ( log ( fabs ( piv_v [ j ] ) / piv_abs ) - markowitz_c * log ( double ( NR [ j ] ) / double ( N_max ) ) ) ;
else
markovitz = fabs ( piv_v [ j ] ) / piv_abs ;
}
else
markovitz = 0 ;
}
else
markovitz = fabs ( piv_v [ j ] ) / piv_abs ;
if ( markovitz > markovitz_max )
{
piv = piv_v [ j ] ;
pivj = pivj_v [ j ] ; //Line number
pivk = pivk_v [ j ] ; //positi
markovitz_max = markovitz ;
}
}
}
else
{
for ( int j = 0 ; j < l ; j + + )
{
if ( N_max > 0 & & NR [ j ] > 0 )
{
if ( fabs ( piv_v [ j ] ) > 0 )
{
if ( markowitz_c > 0 )
markovitz = exp ( log ( fabs ( piv_v [ j ] ) / piv_abs ) - markowitz_c * log ( double ( NR [ j ] ) / double ( N_max ) ) ) ;
else
markovitz = fabs ( piv_v [ j ] ) / piv_abs ;
}
else
markovitz = 0 ;
}
else
markovitz = fabs ( piv_v [ j ] ) / piv_abs ;
if ( NR [ j ] = = 1 )
{
piv = piv_v [ j ] ;
pivj = pivj_v [ j ] ; //Line number
pivk = pivk_v [ j ] ; //positi
markovitz_max = markovitz ;
}
}
}
pivot [ i ] = pivj ;
pivotk [ i ] = pivk ;
pivotv [ i ] = piv ;
line_done [ pivj ] = true ;
/*divide all the non zeros elements of the line pivj by the max_pivot*/
int nb_var = At_Row ( pivj , & first ) ;
for ( int j = 0 ; j < nb_var ; j + + )
{
u [ first - > u_index ] / = piv ;
first = first - > NZE_R_N ;
}
u [ b [ pivj ] ] / = piv ;
/*substract the elements on the non treated lines*/
nb_eq = At_Col ( i , & first ) ;
NonZeroElem * first_piva ;
int nb_var_piva = At_Row ( pivj , & first_piva ) ;
int nb_eq_todo = 0 ;
for ( int j = 0 ; j < nb_eq & & first ; j + + )
{
if ( ! line_done [ first - > r_index ] )
bc [ nb_eq_todo + + ] = first ;
first = first - > NZE_C_N ;
}
//#pragma omp parallel for num_threads(atoi(getenv("DYNARE_NUM_THREADS")))
for ( int j = 0 ; j < nb_eq_todo ; j + + )
{
first = bc [ j ] ;
int row = first - > r_index ;
double first_elem = u [ first - > u_index ] ;
int nb_var_piv = nb_var_piva ;
NonZeroElem * first_piv = first_piva ;
NonZeroElem * first_sub ;
int nb_var_sub = At_Row ( row , & first_sub ) ;
int l_sub = 0 , l_piv = 0 ;
int sub_c_index = first_sub - > c_index , piv_c_index = first_piv - > c_index ;
while ( l_sub < nb_var_sub | | l_piv < nb_var_piv )
{
if ( l_sub < nb_var_sub & & ( sub_c_index < piv_c_index | | l_piv > = nb_var_piv ) )
{
first_sub = first_sub - > NZE_R_N ;
if ( first_sub )
sub_c_index = first_sub - > c_index ;
else
sub_c_index = Size ;
l_sub + + ;
}
else if ( sub_c_index > piv_c_index | | l_sub > = nb_var_sub )
{
int tmp_u_count = Get_u ( ) ;
Insert ( row , first_piv - > c_index , tmp_u_count , 0 ) ;
u [ tmp_u_count ] = - u [ first_piv - > u_index ] * first_elem ;
first_piv = first_piv - > NZE_R_N ;
if ( first_piv )
piv_c_index = first_piv - > c_index ;
else
piv_c_index = Size ;
l_piv + + ;
}
else
{
if ( i = = sub_c_index )
{
firsta = first ;
first_suba = first_sub - > NZE_R_N ;
Delete ( first_sub - > r_index , first_sub - > c_index ) ;
first = firsta - > NZE_C_N ;
first_sub = first_suba ;
if ( first_sub )
sub_c_index = first_sub - > c_index ;
else
sub_c_index = Size ;
l_sub + + ;
first_piv = first_piv - > NZE_R_N ;
if ( first_piv )
piv_c_index = first_piv - > c_index ;
else
piv_c_index = Size ;
l_piv + + ;
}
else
{
u [ first_sub - > u_index ] - = u [ first_piv - > u_index ] * first_elem ;
first_sub = first_sub - > NZE_R_N ;
if ( first_sub )
sub_c_index = first_sub - > c_index ;
else
sub_c_index = Size ;
l_sub + + ;
first_piv = first_piv - > NZE_R_N ;
if ( first_piv )
piv_c_index = first_piv - > c_index ;
else
piv_c_index = Size ;
l_piv + + ;
}
}
}
u [ b [ row ] ] - = u [ b [ pivj ] ] * first_elem ;
}
}
2012-02-17 10:59:39 +01:00
double slowc_lbx = slowc ;
2010-09-24 12:52:58 +02:00
for ( int i = 0 ; i < y_size ; i + + )
ya [ i + it_ * y_size ] = y [ i + it_ * y_size ] ;
slowc_save = slowc ;
2012-02-17 10:59:39 +01:00
simple_bksub ( it_ , Size , slowc_lbx ) ;
2010-09-24 12:52:58 +02:00
End_GE ( Size ) ;
mxFree ( piv_v ) ;
mxFree ( pivj_v ) ;
mxFree ( pivk_v ) ;
mxFree ( NR ) ;
mxFree ( bc ) ;
2012-03-06 11:03:13 +01:00
return false ;
2010-09-24 12:52:58 +02:00
}
void
SparseMatrix : : Solve_ByteCode_Symbolic_Sparse_GaussianElimination ( int Size , bool symbolic , int Block_number )
{
/*Triangularisation at each period of a block using a simple gaussian Elimination*/
t_save_op_s * save_op_s ;
int * save_op = NULL , * save_opa = NULL , * save_opaa = NULL ;
long int nop = 0 , nopa = 0 ;
bool record = false ;
double * piv_v ;
double piv_abs ;
int * pivj_v , * pivk_v , * NR ;
int pivj = 0 , pivk = 0 ;
NonZeroElem * first ;
int tmp_u_count , lag ;
int tbreak = 0 , last_period = periods ;
piv_v = ( double * ) mxMalloc ( Size * sizeof ( double ) ) ;
pivj_v = ( int * ) mxMalloc ( Size * sizeof ( int ) ) ;
pivk_v = ( int * ) mxMalloc ( Size * sizeof ( int ) ) ;
NR = ( int * ) mxMalloc ( Size * sizeof ( int ) ) ;
for ( int t = 0 ; t < periods ; t + + )
{
if ( record & & symbolic )
{
if ( save_op )
{
mxFree ( save_op ) ;
save_op = NULL ;
}
save_op = ( int * ) mxMalloc ( nop * sizeof ( int ) ) ;
nopa = nop ;
}
nop = 0 ;
Clear_u ( ) ;
int ti = t * Size ;
for ( int i = ti ; i < Size + ti ; i + + )
{
/*finding the max-pivot*/
double piv = piv_abs = 0 ;
int nb_eq = At_Col ( i , 0 , & first ) ;
if ( ( symbolic & & t < = start_compare ) | | ! symbolic )
{
int l = 0 , N_max = 0 ;
bool one = false ;
piv_abs = 0 ;
for ( int j = 0 ; j < nb_eq ; j + + )
{
if ( ! line_done [ first - > r_index ] )
{
int k = first - > u_index ;
int jj = first - > r_index ;
int NRow_jj = NRow ( jj ) ;
piv_v [ l ] = u [ k ] ;
double piv_fabs = fabs ( u [ k ] ) ;
pivj_v [ l ] = jj ;
pivk_v [ l ] = k ;
NR [ l ] = NRow_jj ;
if ( NRow_jj = = 1 & & ! one )
{
one = true ;
piv_abs = piv_fabs ;
N_max = NRow_jj ;
}
if ( ! one )
{
if ( piv_fabs > piv_abs )
piv_abs = piv_fabs ;
if ( NRow_jj > N_max )
N_max = NRow_jj ;
}
else
{
if ( NRow_jj = = 1 )
{
if ( piv_fabs > piv_abs )
piv_abs = piv_fabs ;
if ( NRow_jj > N_max )
N_max = NRow_jj ;
}
}
l + + ;
}
first = first - > NZE_C_N ;
}
double markovitz = 0 , markovitz_max = - 9e70 ;
int NR_max = 0 ;
if ( ! one )
{
for ( int j = 0 ; j < l ; j + + )
{
if ( N_max > 0 & & NR [ j ] > 0 )
{
if ( fabs ( piv_v [ j ] ) > 0 )
{
if ( markowitz_c > 0 )
markovitz = exp ( log ( fabs ( piv_v [ j ] ) / piv_abs ) - markowitz_c * log ( double ( NR [ j ] ) / double ( N_max ) ) ) ;
else
markovitz = fabs ( piv_v [ j ] ) / piv_abs ;
}
else
markovitz = 0 ;
}
else
markovitz = fabs ( piv_v [ j ] ) / piv_abs ;
if ( markovitz > markovitz_max )
{
piv = piv_v [ j ] ;
pivj = pivj_v [ j ] ; //Line number
pivk = pivk_v [ j ] ; //positi
markovitz_max = markovitz ;
NR_max = NR [ j ] ;
}
}
}
else
{
for ( int j = 0 ; j < l ; j + + )
{
if ( N_max > 0 & & NR [ j ] > 0 )
{
if ( fabs ( piv_v [ j ] ) > 0 )
{
if ( markowitz_c > 0 )
markovitz = exp ( log ( fabs ( piv_v [ j ] ) / piv_abs ) - markowitz_c * log ( double ( NR [ j ] ) / double ( N_max ) ) ) ;
else
markovitz = fabs ( piv_v [ j ] ) / piv_abs ;
}
else
markovitz = 0 ;
}
else
markovitz = fabs ( piv_v [ j ] ) / piv_abs ;
if ( NR [ j ] = = 1 )
{
piv = piv_v [ j ] ;
pivj = pivj_v [ j ] ; //Line number
pivk = pivk_v [ j ] ; //positi
markovitz_max = markovitz ;
NR_max = NR [ j ] ;
}
}
}
if ( fabs ( piv ) < eps )
2011-02-04 16:53:12 +01:00
mexPrintf ( " ==> Error NR_max=%d, N_max=%d and piv=%f, piv_abs=%f, markovitz_max=%f \n " , NR_max , N_max , piv , piv_abs , markovitz_max ) ;
2010-09-24 12:52:58 +02:00
if ( NR_max = = 0 )
2011-02-04 16:53:12 +01:00
mexPrintf ( " ==> Error NR_max=0 and piv=%f, markovitz_max=%f \n " , piv , markovitz_max ) ;
2010-09-24 12:52:58 +02:00
pivot [ i ] = pivj ;
pivot_save [ i ] = pivj ;
pivotk [ i ] = pivk ;
pivotv [ i ] = piv ;
}
else
{
pivj = pivot [ i - Size ] + Size ;
pivot [ i ] = pivj ;
At_Pos ( pivj , i , & first ) ;
pivk = first - > u_index ;
piv = u [ pivk ] ;
piv_abs = fabs ( piv ) ;
}
line_done [ pivj ] = true ;
if ( symbolic )
{
if ( record )
{
if ( nop + 1 > = nopa )
{
2011-02-04 16:53:12 +01:00
nopa = long ( mem_increasing_factor * ( double ) nopa ) ;
2010-09-24 12:52:58 +02:00
save_op = ( int * ) mxRealloc ( save_op , nopa * sizeof ( int ) ) ;
}
2011-02-04 16:53:12 +01:00
save_op_s = ( t_save_op_s * ) ( & ( save_op [ nop ] ) ) ;
2010-09-24 12:52:58 +02:00
save_op_s - > operat = IFLD ;
save_op_s - > first = pivk ;
save_op_s - > lag = 0 ;
}
nop + = 2 ;
}
if ( piv_abs < eps )
{
ostringstream tmp ;
2011-02-04 16:53:12 +01:00
if ( Block_number > 1 )
2010-09-24 12:52:58 +02:00
tmp < < " in Solve_ByteCode_Symbolic_Sparse_GaussianElimination, singular system in block " < < Block_number + 1 < < " \n " ;
else
2011-02-04 16:53:12 +01:00
tmp < < " in Solve_ByteCode_Symbolic_Sparse_GaussianElimination, singular system \n " ;
2010-09-24 12:52:58 +02:00
throw FatalExceptionHandling ( tmp . str ( ) ) ;
}
/*divide all the non zeros elements of the line pivj by the max_pivot*/
int nb_var = At_Row ( pivj , & first ) ;
NonZeroElem * * bb ;
bb = ( NonZeroElem * * ) mxMalloc ( nb_var * sizeof ( first ) ) ;
for ( int j = 0 ; j < nb_var ; j + + )
{
bb [ j ] = first ;
first = first - > NZE_R_N ;
}
for ( int j = 0 ; j < nb_var ; j + + )
{
first = bb [ j ] ;
u [ first - > u_index ] / = piv ;
if ( symbolic )
{
if ( record )
{
if ( nop + j * 2 + 1 > = nopa )
{
2011-02-04 16:53:12 +01:00
nopa = long ( mem_increasing_factor * ( double ) nopa ) ;
2010-09-24 12:52:58 +02:00
save_op = ( int * ) mxRealloc ( save_op , nopa * sizeof ( int ) ) ;
}
2011-02-04 16:53:12 +01:00
save_op_s = ( t_save_op_s * ) ( & ( save_op [ nop + j * 2 ] ) ) ;
2010-09-24 12:52:58 +02:00
save_op_s - > operat = IFDIV ;
save_op_s - > first = first - > u_index ;
save_op_s - > lag = first - > lag_index ;
}
}
}
mxFree ( bb ) ;
nop + = nb_var * 2 ;
u [ b [ pivj ] ] / = piv ;
if ( symbolic )
{
if ( record )
{
if ( nop + 1 > = nopa )
{
2011-02-04 16:53:12 +01:00
nopa = long ( mem_increasing_factor * ( double ) nopa ) ;
2010-09-24 12:52:58 +02:00
save_op = ( int * ) mxRealloc ( save_op , nopa * sizeof ( int ) ) ;
}
2011-02-04 16:53:12 +01:00
save_op_s = ( t_save_op_s * ) ( & ( save_op [ nop ] ) ) ;
2010-09-24 12:52:58 +02:00
save_op_s - > operat = IFDIV ;
save_op_s - > first = b [ pivj ] ;
save_op_s - > lag = 0 ;
}
nop + = 2 ;
}
/*substract the elements on the non treated lines*/
nb_eq = At_Col ( i , & first ) ;
NonZeroElem * first_piva ;
int nb_var_piva = At_Row ( pivj , & first_piva ) ;
NonZeroElem * * bc ;
bc = ( NonZeroElem * * ) mxMalloc ( nb_eq * sizeof ( first ) ) ;
int nb_eq_todo = 0 ;
for ( int j = 0 ; j < nb_eq & & first ; j + + )
{
if ( ! line_done [ first - > r_index ] )
bc [ nb_eq_todo + + ] = first ;
first = first - > NZE_C_N ;
}
//#pragma omp parallel for num_threads(2) shared(nb_var_piva, first_piva, nopa, nop, save_op, record)
for ( int j = 0 ; j < nb_eq_todo ; j + + )
{
t_save_op_s * save_op_s_l ;
first = bc [ j ] ;
int row = first - > r_index ;
double first_elem = u [ first - > u_index ] ;
if ( symbolic )
{
if ( record )
{
if ( nop + 1 > = nopa )
{
2011-02-04 16:53:12 +01:00
nopa = long ( mem_increasing_factor * ( double ) nopa ) ;
2010-09-24 12:52:58 +02:00
save_op = ( int * ) mxRealloc ( save_op , nopa * sizeof ( int ) ) ;
}
2011-02-04 16:53:12 +01:00
save_op_s_l = ( t_save_op_s * ) ( & ( save_op [ nop ] ) ) ;
2010-09-24 12:52:58 +02:00
save_op_s_l - > operat = IFLD ;
save_op_s_l - > first = first - > u_index ;
save_op_s_l - > lag = abs ( first - > lag_index ) ;
}
nop + = 2 ;
}
int nb_var_piv = nb_var_piva ;
NonZeroElem * first_piv = first_piva ;
NonZeroElem * first_sub ;
int nb_var_sub = At_Row ( row , & first_sub ) ;
int l_sub = 0 ;
int l_piv = 0 ;
int sub_c_index = first_sub - > c_index ;
int piv_c_index = first_piv - > c_index ;
int tmp_lag = first_sub - > lag_index ;
while ( l_sub < nb_var_sub | | l_piv < nb_var_piv )
{
if ( l_sub < nb_var_sub & & ( sub_c_index < piv_c_index | | l_piv > = nb_var_piv ) )
{
//There is no nonzero element at row pivot for this column=> Nothing to do for the current element got to next column
first_sub = first_sub - > NZE_R_N ;
if ( first_sub )
sub_c_index = first_sub - > c_index ;
else
sub_c_index = Size * periods ;
l_sub + + ;
}
else if ( sub_c_index > piv_c_index | | l_sub > = nb_var_sub )
{
// There is an nonzero element at row pivot but not at the current row=> insert a negative element in the current row
tmp_u_count = Get_u ( ) ;
lag = first_piv - > c_index / Size - row / Size ;
//#pragma omp critical
{
Insert ( row , first_piv - > c_index , tmp_u_count , lag ) ;
}
u [ tmp_u_count ] = - u [ first_piv - > u_index ] * first_elem ;
if ( symbolic )
{
if ( record )
{
if ( nop + 2 > = nopa )
{
2011-02-04 16:53:12 +01:00
nopa = long ( mem_increasing_factor * ( double ) nopa ) ;
2010-09-24 12:52:58 +02:00
save_op = ( int * ) mxRealloc ( save_op , nopa * sizeof ( int ) ) ;
}
2011-02-04 16:53:12 +01:00
save_op_s_l = ( t_save_op_s * ) ( & ( save_op [ nop ] ) ) ;
2010-09-24 12:52:58 +02:00
save_op_s_l - > operat = IFLESS ;
save_op_s_l - > first = tmp_u_count ;
save_op_s_l - > second = first_piv - > u_index ;
save_op_s_l - > lag = max ( first_piv - > lag_index , abs ( tmp_lag ) ) ;
}
nop + = 3 ;
}
first_piv = first_piv - > NZE_R_N ;
if ( first_piv )
piv_c_index = first_piv - > c_index ;
else
piv_c_index = Size * periods ;
l_piv + + ;
}
else /*first_sub->c_index==first_piv->c_index*/
{
if ( i = = sub_c_index )
{
NonZeroElem * firsta = first ;
NonZeroElem * first_suba = first_sub - > NZE_R_N ;
Delete ( first_sub - > r_index , first_sub - > c_index ) ;
first = firsta - > NZE_C_N ;
first_sub = first_suba ;
if ( first_sub )
sub_c_index = first_sub - > c_index ;
else
sub_c_index = Size * periods ;
l_sub + + ;
first_piv = first_piv - > NZE_R_N ;
if ( first_piv )
piv_c_index = first_piv - > c_index ;
else
piv_c_index = Size * periods ;
l_piv + + ;
}
else
{
u [ first_sub - > u_index ] - = u [ first_piv - > u_index ] * first_elem ;
if ( symbolic )
{
if ( record )
{
if ( nop + 3 > = nopa )
{
2011-02-04 16:53:12 +01:00
nopa = long ( mem_increasing_factor * ( double ) nopa ) ;
2010-09-24 12:52:58 +02:00
save_op = ( int * ) mxRealloc ( save_op , nopa * sizeof ( int ) ) ;
}
2011-02-04 16:53:12 +01:00
save_op_s_l = ( t_save_op_s * ) ( & ( save_op [ nop ] ) ) ;
2010-09-24 12:52:58 +02:00
save_op_s_l - > operat = IFSUB ;
save_op_s_l - > first = first_sub - > u_index ;
save_op_s_l - > second = first_piv - > u_index ;
save_op_s_l - > lag = max ( abs ( tmp_lag ) , first_piv - > lag_index ) ;
}
nop + = 3 ;
}
first_sub = first_sub - > NZE_R_N ;
if ( first_sub )
sub_c_index = first_sub - > c_index ;
else
sub_c_index = Size * periods ;
l_sub + + ;
first_piv = first_piv - > NZE_R_N ;
if ( first_piv )
piv_c_index = first_piv - > c_index ;
else
piv_c_index = Size * periods ;
l_piv + + ;
}
}
}
u [ b [ row ] ] - = u [ b [ pivj ] ] * first_elem ;
if ( symbolic )
{
if ( record )
{
if ( nop + 3 > = nopa )
{
2011-02-04 16:53:12 +01:00
nopa = long ( mem_increasing_factor * ( double ) nopa ) ;
2010-09-24 12:52:58 +02:00
save_op = ( int * ) mxRealloc ( save_op , nopa * sizeof ( int ) ) ;
}
2011-02-04 16:53:12 +01:00
save_op_s_l = ( t_save_op_s * ) ( & ( save_op [ nop ] ) ) ;
2010-09-24 12:52:58 +02:00
save_op_s_l - > operat = IFSUB ;
save_op_s_l - > first = b [ row ] ;
save_op_s_l - > second = b [ pivj ] ;
save_op_s_l - > lag = abs ( tmp_lag ) ;
}
nop + = 3 ;
}
}
mxFree ( bc ) ;
}
if ( symbolic )
{
if ( record & & ( nop = = nop1 ) )
{
if ( save_opa & & save_opaa )
{
if ( compare ( save_op , save_opa , save_opaa , t , periods , nop , Size ) )
{
tbreak = t ;
tbreak_g = tbreak ;
break ;
}
}
if ( save_opa )
{
if ( save_opaa )
{
mxFree ( save_opaa ) ;
save_opaa = NULL ;
}
save_opaa = ( int * ) mxMalloc ( nop1 * sizeof ( int ) ) ;
memcpy ( save_opaa , save_opa , nop1 * sizeof ( int ) ) ;
}
if ( save_opa )
{
mxFree ( save_opa ) ;
save_opa = NULL ;
}
save_opa = ( int * ) mxMalloc ( nop * sizeof ( int ) ) ;
memcpy ( save_opa , save_op , nop * sizeof ( int ) ) ;
}
else
{
if ( nop = = nop1 )
record = true ;
else
{
record = false ;
if ( save_opa )
{
mxFree ( save_opa ) ;
save_opa = NULL ;
}
if ( save_opaa )
{
mxFree ( save_opaa ) ;
save_opaa = NULL ;
}
}
}
nop2 = nop1 ;
nop1 = nop ;
}
}
mxFree ( piv_v ) ;
mxFree ( pivj_v ) ;
mxFree ( pivk_v ) ;
mxFree ( NR ) ;
nop_all + = nop ;
if ( symbolic )
{
if ( save_op )
mxFree ( save_op ) ;
if ( save_opa )
mxFree ( save_opa ) ;
if ( save_opaa )
mxFree ( save_opaa ) ;
}
/*The backward substitution*/
2012-02-17 10:59:39 +01:00
double slowc_lbx = slowc ;
2010-09-24 12:52:58 +02:00
for ( int i = 0 ; i < y_size * ( periods + y_kmin ) ; i + + )
ya [ i ] = y [ i ] ;
slowc_save = slowc ;
2012-02-17 10:59:39 +01:00
bksub ( tbreak , last_period , Size , slowc_lbx ) ;
2010-09-24 12:52:58 +02:00
End_GE ( Size ) ;
}
2012-03-06 11:03:13 +01:00
bool
2010-10-18 17:28:21 +02:00
SparseMatrix : : Simulate_Newton_One_Boundary ( 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 stack_solve_algo , int solve_algo )
2010-09-24 12:52:58 +02:00
{
int i , j ;
2011-01-03 14:15:12 +01:00
mxArray * b_m = NULL , * A_m = NULL , * x0_m = NULL ;
2010-09-24 12:52:58 +02:00
Clear_u ( ) ;
error_not_printed = true ;
2012-03-06 11:03:13 +01:00
bool singular_system = false ;
2010-09-24 12:52:58 +02:00
u_count_alloc_save = u_count_alloc ;
2011-02-04 16:53:12 +01:00
if ( isnan ( res1 ) | | isinf ( res1 ) | | ( res2 > 12 * g0 & & iter > 0 ) )
2010-09-24 12:52:58 +02:00
{
if ( iter = = 0 | | fabs ( slowc_save ) < 1e-8 )
{
for ( j = 0 ; j < y_size ; j + + )
{
2012-06-08 19:00:57 +02:00
# ifdef DEBUG
2010-09-24 12:52:58 +02:00
bool select = false ;
2012-06-08 19:00:57 +02:00
# endif
2010-09-24 12:52:58 +02:00
for ( int i = 0 ; i < Size ; i + + )
if ( j = = index_vara [ i ] )
{
2012-06-08 19:00:57 +02:00
# ifdef DEBUG
2010-09-24 12:52:58 +02:00
select = true ;
2012-06-08 19:00:57 +02:00
# endif
2010-09-24 12:52:58 +02:00
break ;
}
2012-05-20 14:37:30 +02:00
# ifdef DEBUG
2010-09-24 12:52:58 +02:00
if ( select )
2010-11-20 15:45:15 +01:00
mexPrintf ( " -> variable %s (%d) at time %d = %f direction = %f \n " , get_variable ( eEndogenous , j ) . c_str ( ) , j + 1 , it_ , y [ j + it_ * y_size ] , direction [ j + it_ * y_size ] ) ;
2010-09-24 12:52:58 +02:00
else
2010-11-20 15:45:15 +01:00
mexPrintf ( " variable %s (%d) at time %d = %f direction = %f \n " , get_variable ( eEndogenous , j ) . c_str ( ) , j + 1 , it_ , y [ j + it_ * y_size ] , direction [ j + it_ * y_size ] ) ;
2012-05-20 14:37:30 +02:00
# endif
2010-09-24 12:52:58 +02:00
}
if ( steady_state )
{
if ( iter = = 0 )
mexPrintf ( " the initial values of endogenous variables are too far from the solution. \n Change them! \n " ) ;
else
mexPrintf ( " dynare cannot improve the simulation in block %d at time %d (variable %d) \n " , blck + 1 , it_ + 1 , index_vara [ max_res_idx ] + 1 ) ;
mexEvalString ( " drawnow; " ) ;
2012-03-06 11:03:13 +01:00
return singular_system ;
2010-09-24 12:52:58 +02:00
}
else
{
ostringstream tmp ;
if ( iter = = 0 )
tmp < < " in Simulate_Newton_One_Boundary, The initial values of endogenous variables are too far from the solution. \n Change them! \n " ;
else
tmp < < " in Simulate_Newton_One_Boundary, Dynare cannot improve the simulation in block " < < blck + 1 < < " at time " < < it_ + 1 < < " (variable " < < index_vara [ max_res_idx ] + 1 < < " %d) \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
}
}
2011-02-04 16:53:12 +01:00
if ( ! ( isnan ( res1 ) | | isinf ( res1 ) ) & & ! ( isnan ( g0 ) | | isinf ( g0 ) ) )
2010-09-24 12:52:58 +02:00
{
if ( try_at_iteration = = 0 )
{
prev_slowc_save = slowc_save ;
2011-02-04 16:53:12 +01:00
slowc_save = max ( - gp0 / ( 2 * ( res2 - g0 - gp0 ) ) , 0.1 ) ;
2010-09-24 12:52:58 +02:00
}
else
{
double t1 = res2 - gp0 * slowc_save - g0 ;
double t2 = glambda2 - gp0 * prev_slowc_save - g0 ;
double a = ( 1 / ( slowc_save * slowc_save ) * t1 - 1 / ( prev_slowc_save * prev_slowc_save ) * t2 ) / ( slowc_save - prev_slowc_save ) ;
double b = ( - prev_slowc_save / ( slowc_save * slowc_save ) * t1 + slowc_save / ( prev_slowc_save * prev_slowc_save ) * t2 ) / ( slowc_save - prev_slowc_save ) ;
prev_slowc_save = slowc_save ;
2011-02-04 16:53:12 +01:00
slowc_save = max ( min ( - b + sqrt ( b * b - 3 * a * gp0 ) / ( 3 * a ) , 0.5 * slowc_save ) , 0.1 * slowc_save ) ;
2010-09-24 12:52:58 +02:00
}
glambda2 = res2 ;
2011-02-04 16:53:12 +01:00
try_at_iteration + + ;
2010-09-24 12:52:58 +02:00
}
else
{
prev_slowc_save = slowc_save ;
slowc_save / = 1.1 ;
}
if ( print_it )
mexPrintf ( " Error: Simulation diverging, trying to correct it using slowc=%f \n " , slowc_save ) ;
for ( i = 0 ; i < y_size ; i + + )
y [ i + it_ * y_size ] = ya [ i + it_ * y_size ] + slowc_save * direction [ i + it_ * y_size ] ;
iter - - ;
2012-03-06 11:03:13 +01:00
return singular_system ;
2010-09-24 12:52:58 +02:00
}
if ( cvg )
{
2012-03-06 11:03:13 +01:00
return singular_system ;
2010-09-24 12:52:58 +02:00
}
2011-02-04 16:53:12 +01:00
if ( print_it )
2010-09-24 12:52:58 +02:00
{
2010-10-18 17:28:21 +02:00
//mexPrintf("solwc=%f g0=%f res2=%f glambda2=%f\n",slowc_save,g0, res2, glambda2);
if ( steady_state )
{
switch ( solve_algo )
{
2011-02-04 16:53:12 +01:00
case 0 :
mexPrintf ( " MODEL STEADY STATE: MATLAB fsolve \n " ) ;
break ;
case 1 :
mexPrintf ( " MODEL STEADY STATE: MATLAB solve1 \n " ) ;
break ;
case 2 :
case 4 :
mexPrintf ( " MODEL STEADY STATE: block decomposition + MATLAB solve1 \n " ) ;
break ;
case 3 :
mexPrintf ( " MODEL STEADY STATE: MATLAB csolve \n " ) ;
break ;
case 5 :
mexPrintf ( " MODEL STEADY STATE: (method=ByteCode own solver) \n " ) ;
break ;
case 6 :
mexPrintf ( " MODEL STEADY STATE: Sparse LU \n " ) ;
break ;
case 7 :
mexPrintf ( " MODEL STEADY STATE: (method=GMRES) \n " ) ;
break ;
case 8 :
mexPrintf ( " MODEL STEADY STATE: (method=BiCGStab) \n " ) ;
break ;
default :
mexPrintf ( " MODEL STEADY STATE: (method=Unknown - %d - ) \n " , stack_solve_algo ) ;
2010-10-18 17:28:21 +02:00
}
}
2010-09-24 12:52:58 +02:00
mexPrintf ( " ----------------------------------- \n " ) ;
mexPrintf ( " Simulate iteration no %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 " ) ;
}
2010-10-18 17:28:21 +02:00
bool zero_solution ;
2010-10-22 16:49:47 +02:00
if ( ( solve_algo = = 5 & & steady_state ) | | ( stack_solve_algo = = 5 & & ! steady_state ) )
2012-03-06 11:03:13 +01:00
Simple_Init ( Size , IM_i , zero_solution ) ;
2010-09-24 12:52:58 +02:00
else
{
2011-02-04 16:53:12 +01:00
b_m = mxCreateDoubleMatrix ( Size , 1 , mxREAL ) ;
2010-09-24 12:52:58 +02:00
if ( ! b_m )
{
ostringstream tmp ;
tmp < < " in Simulate_Newton_One_Boundary, can't allocate b_m vector \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
}
2011-02-04 16:53:12 +01:00
A_m = mxCreateSparse ( Size , Size , min ( int ( IM_i . size ( ) * 2 ) , Size * Size ) , mxREAL ) ;
2010-09-24 12:52:58 +02:00
if ( ! A_m )
{
ostringstream tmp ;
tmp < < " in Simulate_Newton_One_Boundary, can't allocate A_m matrix \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
}
2011-02-04 16:53:12 +01:00
x0_m = mxCreateDoubleMatrix ( Size , 1 , mxREAL ) ;
2011-01-03 14:15:12 +01:00
if ( ! x0_m )
{
ostringstream tmp ;
tmp < < " in Simulate_Newton_One_Boundary, can't allocate x0_m vector \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
}
Init_Matlab_Sparse_Simple ( Size , IM_i , A_m , b_m , zero_solution , x0_m ) ;
2010-10-18 17:28:21 +02:00
}
if ( zero_solution )
{
for ( int i = 0 ; i < Size ; i + + )
{
int eq = index_vara [ i ] ;
2011-02-04 16:53:12 +01:00
double yy = - ( y [ eq + it_ * y_size ] ) ;
2010-10-18 17:28:21 +02:00
direction [ eq ] = yy ;
y [ eq + it_ * y_size ] + = slowc * yy ;
}
}
else
{
2010-10-22 16:49:47 +02:00
if ( ( solve_algo = = 5 & & steady_state ) | | ( stack_solve_algo = = 5 & & ! steady_state ) )
2012-03-06 11:03:13 +01:00
singular_system = Solve_ByteCode_Sparse_GaussianElimination ( Size , blck , steady_state , it_ ) ;
2010-10-22 16:49:47 +02:00
else if ( ( solve_algo = = 7 & & steady_state ) | | ( stack_solve_algo = = 2 & & ! steady_state ) )
2011-01-03 14:15:12 +01:00
Solve_Matlab_GMRES ( A_m , b_m , Size , slowc , blck , false , it_ , steady_state , x0_m ) ;
2010-10-22 16:49:47 +02:00
else if ( ( solve_algo = = 8 & & steady_state ) | | ( stack_solve_algo = = 3 & & ! steady_state ) )
2011-01-03 14:15:12 +01:00
Solve_Matlab_BiCGStab ( A_m , b_m , Size , slowc , blck , false , it_ , x0_m , steady_state ) ;
2010-10-22 16:49:47 +02:00
else if ( ( solve_algo = = 6 & & steady_state ) | | ( ( stack_solve_algo = = 0 | | stack_solve_algo = = 1 ) & & ! steady_state ) )
2010-10-18 17:28:21 +02:00
Solve_Matlab_LU_UMFPack ( A_m , b_m , Size , slowc , false , it_ ) ;
2010-09-24 12:52:58 +02:00
}
2012-03-06 11:03:13 +01:00
return singular_system ;
2010-09-24 12:52:58 +02:00
}
void
SparseMatrix : : Simulate_Newton_Two_Boundaries ( 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 stack_solve_algo , unsigned int endo_name_length , char * P_endo_names )
{
2009-08-25 11:43:01 +02:00
if ( start_compare = = 0 )
start_compare = y_kmin ;
u_count_alloc_save = u_count_alloc ;
2007-11-21 00:27:30 +01:00
clock_t t1 = clock ( ) ;
2009-07-10 17:10:11 +02:00
nop1 = 0 ;
2009-08-25 11:43:01 +02:00
error_not_printed = true ;
2011-01-03 14:15:12 +01:00
mxArray * b_m = NULL , * A_m = NULL , * x0_m = NULL ;
2009-08-25 11:43:01 +02:00
if ( iter > 0 )
2008-12-19 11:24:31 +01:00
{
2011-12-12 12:40:46 +01:00
if ( print_it )
{
mexPrintf ( " Sim : %f ms \n " , ( 1000.0 * ( double ( clock ( ) ) - double ( time00 ) ) ) / double ( CLOCKS_PER_SEC ) ) ;
mexEvalString ( " drawnow; " ) ;
}
2009-08-25 11:43:01 +02:00
time00 = clock ( ) ;
2008-12-19 11:24:31 +01:00
}
2011-02-04 16:53:12 +01:00
if ( isnan ( res1 ) | | isinf ( res1 ) | | ( res2 > 12 * g0 & & iter > 0 ) )
2010-02-06 15:07:56 +01:00
{
2010-09-24 12:52:58 +02:00
if ( iter = = 0 | | fabs ( slowc_save ) < 1e-8 )
2010-02-06 15:07:56 +01:00
{
2011-02-04 16:53:12 +01:00
for ( int j = 0 ; j < y_size ; j + + )
2010-02-06 15:07:56 +01:00
{
2010-08-18 13:51:57 +02:00
ostringstream res ;
for ( unsigned int i = 0 ; i < endo_name_length ; i + + )
2010-10-18 17:28:21 +02:00
if ( P_endo_names [ CHAR_LENGTH * ( j + i * y_size ) ] ! = ' ' )
res < < P_endo_names [ CHAR_LENGTH * ( j + i * y_size ) ] ;
2010-02-06 15:07:56 +01:00
bool select = false ;
for ( int i = 0 ; i < Size ; i + + )
if ( j = = index_vara [ i ] )
{
select = true ;
break ;
}
if ( select )
2010-08-18 13:51:57 +02:00
mexPrintf ( " -> variable %s (%d) at time %d = %f direction = %f \n " , res . str ( ) . c_str ( ) , j + 1 , it_ , y [ j + it_ * y_size ] , direction [ j + it_ * y_size ] ) ;
2010-02-06 15:07:56 +01:00
else
2010-08-18 13:51:57 +02:00
mexPrintf ( " variable %s (%d) at time %d = %f direction = %f \n " , res . str ( ) . c_str ( ) , j + 1 , it_ , y [ j + it_ * y_size ] , direction [ j + it_ * y_size ] ) ;
2010-02-06 15:07:56 +01:00
}
2010-09-24 12:52:58 +02:00
ostringstream Error ;
2011-02-04 16:53:12 +01:00
if ( iter = = 0 )
2010-09-24 12:52:58 +02:00
Error < < " in Simulate_Newton_Two_Boundaries, the initial values of endogenous variables are too far from the solution. \n Change them! \n " ;
else
Error < < " in Simulate_Newton_Two_Boundaries, dynare cannot improve the simulation in block " < < blck + 1 < < " at time " < < it_ + 1 < < " (variable " < < index_vara [ max_res_idx ] + 1 < < " ) \n " ;
//Error << filename << " stopped";
throw FatalExceptionHandling ( Error . str ( ) ) ;
2010-02-06 15:07:56 +01:00
}
2011-02-04 16:53:12 +01:00
if ( ! ( isnan ( res1 ) | | isinf ( res1 ) ) & & ! ( isnan ( g0 ) | | isinf ( g0 ) ) & & ( stack_solve_algo = = 4 | | stack_solve_algo = = 5 ) )
2010-02-06 15:07:56 +01:00
{
if ( try_at_iteration = = 0 )
{
prev_slowc_save = slowc_save ;
2011-02-04 16:53:12 +01:00
slowc_save = max ( - gp0 / ( 2 * ( res2 - g0 - gp0 ) ) , 0.1 ) ;
2010-02-06 15:07:56 +01:00
}
else
{
double t1 = res2 - gp0 * slowc_save - g0 ;
double t2 = glambda2 - gp0 * prev_slowc_save - g0 ;
double a = ( 1 / ( slowc_save * slowc_save ) * t1 - 1 / ( prev_slowc_save * prev_slowc_save ) * t2 ) / ( slowc_save - prev_slowc_save ) ;
double b = ( - prev_slowc_save / ( slowc_save * slowc_save ) * t1 + slowc_save / ( prev_slowc_save * prev_slowc_save ) * t2 ) / ( slowc_save - prev_slowc_save ) ;
prev_slowc_save = slowc_save ;
2011-02-04 16:53:12 +01:00
slowc_save = max ( min ( - b + sqrt ( b * b - 3 * a * gp0 ) / ( 3 * a ) , 0.5 * slowc_save ) , 0.1 * slowc_save ) ;
2010-02-06 15:07:56 +01:00
}
glambda2 = res2 ;
2011-02-04 16:53:12 +01:00
try_at_iteration + + ;
if ( slowc_save < = 0.1 )
2010-07-23 11:20:24 +02:00
{
2010-09-24 12:52:58 +02:00
for ( int i = 0 ; i < y_size * ( periods + y_kmin ) ; i + + )
2010-07-23 11:20:24 +02:00
y [ i ] = ya [ i ] + direction [ i ] ;
g0 = res2 ;
gp0 = - res2 ;
try_at_iteration = 0 ;
iter - - ;
2010-09-24 12:52:58 +02:00
return ;
2010-07-23 11:20:24 +02:00
}
2010-02-06 15:07:56 +01:00
}
else
{
prev_slowc_save = slowc_save ;
slowc_save / = 1.1 ;
}
if ( print_it )
2010-08-18 13:51:57 +02:00
{
if ( isnan ( res1 ) | | isinf ( res1 ) )
2010-09-24 12:52:58 +02:00
mexPrintf ( " The model cannot be evaluated, trying to correct it using slowc=%f \n " , slowc_save ) ;
2010-08-18 13:51:57 +02:00
else
2010-09-24 12:52:58 +02:00
mexPrintf ( " Simulation diverging, trying to correct it using slowc=%f \n " , slowc_save ) ;
2010-08-18 13:51:57 +02:00
}
2010-09-24 12:52:58 +02:00
for ( int i = 0 ; i < y_size * ( periods + y_kmin ) ; i + + )
2009-08-25 11:43:01 +02:00
y [ i ] = ya [ i ] + slowc_save * direction [ i ] ;
2007-11-21 00:27:30 +01:00
iter - - ;
2010-09-24 12:52:58 +02:00
return ;
2007-11-21 00:27:30 +01:00
}
2010-09-24 12:52:58 +02:00
2010-07-23 11:20:24 +02:00
u_count + = u_count_init ;
if ( stack_solve_algo = = 5 )
2007-11-21 00:27:30 +01:00
{
2010-07-23 11:20:24 +02:00
if ( alt_symbolic & & alt_symbolic_count < alt_symbolic_count_max )
{
mexPrintf ( " Pivoting method will be applied only to the first periods. \n " ) ;
alt_symbolic = false ;
symbolic = true ;
markowitz_c = markowitz_c_s ;
alt_symbolic_count + + ;
}
if ( ( ( res1 / res1a - 1 ) > - 0.3 ) & & symbolic & & iter > 0 )
2009-08-25 11:43:01 +02:00
{
2010-07-23 11:20:24 +02:00
if ( restart > 2 )
{
mexPrintf ( " Divergence or slowdown occured during simulation. \n In the next iteration, pivoting method will be applied to all periods. \n " ) ;
symbolic = false ;
alt_symbolic = true ;
markowitz_c_s = markowitz_c ;
markowitz_c = 0 ;
}
else
{
mexPrintf ( " Divergence or slowdown occured during simulation. \n In the next iteration, pivoting method will be applied for a longer period. \n " ) ;
start_compare = min ( tbreak_g , periods ) ;
restart + + ;
}
2009-08-25 11:43:01 +02:00
}
2009-07-21 17:50:12 +02:00
else
2007-11-21 00:27:30 +01:00
{
2010-07-23 11:20:24 +02:00
start_compare = max ( y_kmin , minimal_solving_periods ) ;
restart = 0 ;
2007-11-21 00:27:30 +01:00
}
}
2009-08-25 11:43:01 +02:00
res1a = res1 ;
2009-07-21 17:50:12 +02:00
2009-08-25 11:43:01 +02:00
if ( print_it )
2008-12-19 11:24:31 +01:00
{
2010-07-23 11:20:24 +02:00
if ( iter = = 0 )
{
2011-02-04 16:53:12 +01:00
switch ( stack_solve_algo )
{
case 0 :
mexPrintf ( " MODEL SIMULATION: (method=Sparse LU) \n " ) ;
break ;
case 1 :
mexPrintf ( " MODEL SIMULATION: (method=Relaxation) \n " ) ;
break ;
case 2 :
mexPrintf ( " MODEL SIMULATION: (method=GMRES) \n " ) ;
break ;
case 3 :
mexPrintf ( " MODEL SIMULATION: (method=BiCGStab) \n " ) ;
break ;
case 4 :
mexPrintf ( " MODEL SIMULATION: (method=Sparse LU & optimal path length) \n " ) ;
break ;
case 5 :
mexPrintf ( " MODEL SIMULATION: (method=ByteCode own solver) \n " ) ;
break ;
default :
mexPrintf ( " MODEL SIMULATION: (method=Unknown - %d - ) \n " , stack_solve_algo ) ;
2010-07-23 11:20:24 +02:00
}
}
2008-12-19 11:24:31 +01:00
mexPrintf ( " ----------------------------------- \n " ) ;
2009-12-16 14:21:31 +01:00
mexPrintf ( " Simulate iteration no %d \n " , iter + 1 ) ;
2009-08-25 11:43:01 +02:00
mexPrintf ( " max. error=%.10e \n " , double ( max_res ) ) ;
mexPrintf ( " sqr. error=%.10e \n " , double ( res2 ) ) ;
mexPrintf ( " abs. error=%.10e \n " , double ( res1 ) ) ;
2008-12-19 11:24:31 +01:00
mexPrintf ( " ----------------------------------- \n " ) ;
2010-07-23 11:20:24 +02:00
mexEvalString ( " drawnow; " ) ;
2008-12-19 11:24:31 +01:00
}
2007-11-21 00:27:30 +01:00
if ( cvg )
2009-07-10 17:10:11 +02:00
{
2010-09-24 12:52:58 +02:00
return ;
2007-11-21 00:27:30 +01:00
}
else
2009-08-25 11:43:01 +02:00
{
2010-07-23 11:20:24 +02:00
if ( stack_solve_algo = = 5 )
Init_GE ( periods , y_kmin , y_kmax , Size , IM_i ) ;
else
2009-08-25 11:43:01 +02:00
{
2011-02-04 16:53:12 +01:00
b_m = mxCreateDoubleMatrix ( periods * Size , 1 , mxREAL ) ;
2010-07-23 11:20:24 +02:00
if ( ! b_m )
2009-08-25 11:43:01 +02:00
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Simulate_Newton_Two_Boundaries, can't allocate b_m vector \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2009-08-25 11:43:01 +02:00
}
2011-02-04 16:53:12 +01:00
x0_m = mxCreateDoubleMatrix ( periods * Size , 1 , mxREAL ) ;
2011-01-03 14:15:12 +01:00
if ( ! x0_m )
{
ostringstream tmp ;
tmp < < " in Simulate_Newton_Two_Boundaries, can't allocate x0_m vector \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
}
2010-07-23 11:20:24 +02:00
A_m = mxCreateSparse ( periods * Size , periods * Size , IM_i . size ( ) * periods * 2 , mxREAL ) ;
if ( ! A_m )
2009-08-25 11:43:01 +02:00
{
2010-09-24 12:52:58 +02:00
ostringstream tmp ;
tmp < < " in Simulate_Newton_Two_Boundaries, can't allocate A_m matrix \n " ;
throw FatalExceptionHandling ( tmp . str ( ) ) ;
2010-07-23 11:20:24 +02:00
}
2011-01-03 14:15:12 +01:00
Init_Matlab_Sparse ( periods , y_kmin , y_kmax , Size , IM_i , A_m , b_m , x0_m ) ;
2010-07-23 11:20:24 +02:00
}
2010-09-24 12:52:58 +02:00
if ( stack_solve_algo = = 0 | | stack_solve_algo = = 4 )
Solve_Matlab_LU_UMFPack ( A_m , b_m , Size , slowc , true , 0 ) ;
2010-09-17 09:57:38 +02:00
else if ( stack_solve_algo = = 1 )
Solve_Matlab_Relaxation ( A_m , b_m , Size , slowc , true , 0 ) ;
2010-07-23 11:20:24 +02:00
else if ( stack_solve_algo = = 2 )
2011-01-03 14:15:12 +01:00
Solve_Matlab_GMRES ( A_m , b_m , Size , slowc , blck , true , 0 , false , x0_m ) ;
2010-07-23 11:20:24 +02:00
else if ( stack_solve_algo = = 3 )
2011-01-03 14:15:12 +01:00
Solve_Matlab_BiCGStab ( A_m , b_m , Size , slowc , blck , true , 0 , x0_m , false ) ;
2010-09-24 12:52:58 +02:00
else if ( stack_solve_algo = = 5 )
Solve_ByteCode_Symbolic_Sparse_GaussianElimination ( Size , symbolic , blck ) ;
2009-08-25 11:43:01 +02:00
}
2007-11-21 00:27:30 +01:00
if ( print_it )
{
clock_t t2 = clock ( ) ;
2009-08-25 11:43:01 +02:00
mexPrintf ( " (** %f milliseconds **) \n " , 1000.0 * ( double ( t2 ) - double ( t1 ) ) / double ( CLOCKS_PER_SEC ) ) ;
2007-11-21 00:27:30 +01:00
mexEvalString ( " drawnow; " ) ;
}
2009-08-25 11:43:01 +02:00
time00 = clock ( ) ;
if ( tbreak_g = = 0 )
tbreak_g = periods ;
2010-09-24 12:52:58 +02:00
return ;
2007-11-21 00:27:30 +01:00
}
void
SparseMatrix : : fixe_u ( double * * u , int u_count_int , int max_lag_plus_max_lead_plus_1 )
{
2009-08-25 11:43:01 +02:00
u_count = u_count_int * periods ;
2007-11-21 00:27:30 +01:00
u_count_alloc = 2 * u_count ;
2010-12-31 16:29:39 +01:00
# ifdef DEBUG
2011-02-04 16:53:12 +01:00
mexPrintf ( " fixe_u : alloc(%d double) \n " , u_count_alloc ) ;
2010-11-20 15:45:15 +01:00
# endif
2009-08-25 11:43:01 +02:00
( * u ) = ( double * ) mxMalloc ( u_count_alloc * sizeof ( double ) ) ;
2010-12-31 16:29:39 +01:00
# ifdef DEBUG
2011-02-04 16:53:12 +01:00
mexPrintf ( " *u=%d \n " , * u ) ;
2010-11-20 15:45:15 +01:00
# endif
2007-11-21 00:27:30 +01:00
memset ( ( * u ) , 0 , u_count_alloc * sizeof ( double ) ) ;
2009-08-25 11:43:01 +02:00
u_count_init = max_lag_plus_max_lead_plus_1 ;
2007-11-21 00:27:30 +01:00
}