Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • 4.3
  • 4.4
  • 4.5
  • 4.6
  • 5.x
  • 6.x
  • asm
  • aux_func
  • clang+openmp
  • dates-and-dseries-improvements
  • declare_vars_in_model_block
  • dmm
  • dragonfly
  • dynare_minreal
  • eigen
  • error_msg_undeclared_model_vars
  • estim_params
  • exo_steady_state
  • gpm-optimal-policy
  • julia
  • madysson
  • master
  • mex-GetPowerDeriv
  • penalty
  • separateM_
  • slice
  • sphinx-doc-experimental
  • static_aux_vars
  • time-varying-information-set
  • various_fixes
  • 3.062
  • 3.063
  • 4.0.0
  • 4.0.1
  • 4.0.2
  • 4.0.3
  • 4.0.4
  • 4.1-alpha1
  • 4.1-alpha2
  • 4.1.0
  • 4.1.1
  • 4.1.2
  • 4.1.3
  • 4.2.0
  • 4.2.1
  • 4.2.2
  • 4.2.3
  • 4.2.4
  • 4.2.5
  • 4.3.0
  • 4.3.1
  • 4.3.2
  • 4.3.3
  • 4.4-beta1
  • 4.4.0
  • 4.4.1
  • 4.4.2
  • 4.4.3
  • 4.5.0
  • 4.5.1
  • 4.5.2
  • 4.5.3
  • 4.5.4
  • 4.5.5
  • 4.5.6
  • 4.5.7
  • 4.6-beta1
  • 4.6.0
  • 4.6.0-rc1
  • 4.6.0-rc2
  • 4.6.1
  • 4.6.2
  • 4.6.3
  • 4.6.4
  • 4.7-beta1
  • 4.7-beta2
  • 4.7-beta3
  • 5.0
  • 5.0-rc1
  • 5.1
  • 5.2
  • 5.3
  • 5.4
  • 5.5
  • 6-beta1
  • 6-beta2
  • 6.0
  • 6.1
  • 6.2
  • 6.3
90 results

Target

Select target project
  • giovanma/dynare
  • giorgiomas/dynare
  • Vermandel/dynare
  • Dynare/dynare
  • normann/dynare
  • MichelJuillard/dynare
  • wmutschl/dynare
  • FerhatMihoubi/dynare
  • sebastien/dynare
  • lnsongxf/dynare
  • rattoma/dynare
  • CIMERS/dynare
  • FredericKarame/dynare
  • SumuduK/dynare
  • MinjeJeon/dynare
  • camilomrch/dynare
  • DoraK/dynare
  • avtishin/dynare
  • selma/dynare
  • claudio_olguin/dynare
  • jeffjiang07/dynare
  • EthanSystem/dynare
  • stepan-a/dynare
  • wjgatt/dynare
  • JohannesPfeifer/dynare
  • gboehl/dynare
  • chskcau/dynare-doc-fixes
27 results
Select Git revision
  • 4.3
  • 4.4
  • 4.5
  • 4.6
  • 4.7
  • 5.x
  • 6.x
  • Andreasen
  • DSMH
  • aux_func
  • dates-and-dseries-improvements
  • declare_vars_in_model_block
  • dmm
  • doc_bib
  • dsge_Base_test
  • dynamic-striated
  • eigen
  • error_msg_undeclared_model_vars
  • estim_params
  • exceptions
  • exist
  • exo_steady_state
  • filter_initial_state
  • gmhmaxlik
  • gpm-optimal-policy
  • julia
  • manual_error
  • master
  • merge-initvalfile-fix
  • mex
  • mex-GetPowerDeriv
  • new_ep
  • occbin
  • osr_analytic
  • penalty
  • remove-@dates
  • remove-@dseries
  • remove-utilities-tests
  • rmExtraExo
  • semi_doc
  • separateM_
  • slice
  • sphinx-doc-experimental
  • static_aux_vars
  • temporary_terms
  • trust_region_legacy
  • 3.062
  • 3.063
  • 4.0.0
  • 4.0.1
  • 4.0.2
  • 4.0.3
  • 4.0.4
  • 4.1-alpha1
  • 4.1-alpha2
  • 4.1.0
  • 4.1.1
  • 4.1.2
  • 4.1.3
  • 4.2.0
  • 4.2.1
  • 4.2.2
  • 4.2.3
  • 4.2.4
  • 4.2.5
  • 4.3.0
  • 4.3.1
  • 4.3.2
  • 4.3.3
  • 4.4-beta1
  • 4.4.0
  • 4.4.1
  • 4.4.2
  • 4.4.3
  • 4.5.0
  • 4.5.1
  • 4.5.2
  • 4.5.3
  • 4.5.4
  • 4.5.5
  • 4.5.6
81 results
Show changes
Show whitespace changes
Inline Side-by-side
Showing
with 0 additions and 3091 deletions
dynare++/src/dynglob.lex deleted 100644 → 0
View file @ 0691b083
%{
#include "parser/cc/location.h"
#include "dynglob_tab.hh"
extern YYLTYPE dynglob_lloc;
#define YY_USER_ACTION SET_LLOC(dynglob_);
%}
%option nounput
%option noyy_top_state
%option stack
%option yylineno
%option prefix="dynglob_"
%option never-interactive
%x CMT
%%
/* comments */
<*>"/*" {yy_push_state(CMT);}
<CMT>[^*\n]*
<CMT>"*"+[^*/\n]*
<CMT>"*"+"/" {yy_pop_state();}
<CMT>[\n]
"//".*\n
/* initial spaces or tabs are ignored */
[ \t\r\n\0]
var {return VAR;}
varexo {return VAREXO;}
parameters {return PARAMETERS;}
model {return MODEL;}
end {return END;}
initval {return INITVAL;}
order {return ORDER;}
vcov {return VCOV;}
planner_objective {return PLANNEROBJECTIVE;}
planner_discount {return PLANNERDISCOUNT;}
/* names */
[A-Za-z_][A-Za-z0-9_]* {
dynglob_lval.string = dynglob_text;
return NAME;
}
; {return SEMICOLON;}
, {return COMMA;}
= {return EQUAL_SIGN;}
\[ {return LEFT_BRACKET;}
\] {return RIGHT_BRACKET;}
. {
dynglob_lval.character = dynglob_text[0];
return CHARACTER;
}
%%
int dynglob_wrap()
{
return 1;
}
void dynglob__destroy_buffer(void* p)
{
dynglob__delete_buffer((YY_BUFFER_STATE)p);
}
dynare++/src/dynglob.y deleted 100644 → 0
View file @ 0691b083
%{
// Copyright (C) 2006-2011, Ondra Kamenik
#include "parser/cc/location.h"
#include "dynare_model.h"
#include "dynglob_tab.hh"
#include <stdio.h>
void dynglob_error(const char*);
int dynglob_lex(void);
extern int dynglob_lineno;
extern ogdyn::DynareParser* dynare_parser;
int symblist_flag;
// static void print_token_value1 (FILE *, int, YYSTYPE);
//#define YYPRINT(file, type, value) print_token_value1 (file, type, value)
%}
%union {
int integer;
char *string;
char character;
}
%token END INITVAL MODEL PARAMETERS VAR VAREXO SEMICOLON COMMA EQUAL_SIGN CHARACTER
%token VCOV LEFT_BRACKET RIGHT_BRACKET ORDER PLANNEROBJECTIVE PLANNERDISCOUNT
%token <string> NAME;
%name-prefix="dynglob_"
%locations
%error-verbose
%%
dynare_file : preamble paramset model rest {
dynare_parser->set_paramset_pos(@2.off, @3.off);}
| preamble model rest {
dynare_parser->set_paramset_pos(0, 0);}
| preamble paramset planner model rest {
dynare_parser->set_paramset_pos(@2.off, @3.off);}
;
preamble : preamble preamble_statement | preamble_statement;
preamble_statement : var | varexo | parameters;
var : VAR {symblist_flag=1;} symblist SEMICOLON;
varexo : VAREXO {symblist_flag=2;} symblist SEMICOLON;
parameters : PARAMETERS {symblist_flag=3;} symblist SEMICOLON;
symblist : symblist NAME {dynare_parser->add_name($2,symblist_flag);}
| symblist COMMA NAME {dynare_parser->add_name($3,symblist_flag);}
| NAME {dynare_parser->add_name($1,symblist_flag);}
;
paramset : recnameset;
recnameset : recnameset onenameset | onenameset;
onenameset : NAME EQUAL_SIGN material SEMICOLON;
material : material CHARACTER | material NAME | NAME | CHARACTER;
model : MODEL SEMICOLON equations END SEMICOLON {
dynare_parser->set_model_pos(@3.off, @4.off);
};
equations : equations equation | equation;
equation : material EQUAL_SIGN material SEMICOLON | material SEMICOLON;
rest : rest_statement | rest rest_statement;
rest_statement : initval | vcov | order | planner;
initval : INITVAL SEMICOLON recnameset END SEMICOLON {
dynare_parser->set_initval_pos(@3.off, @4.off);
};
vcov : VCOV EQUAL_SIGN LEFT_BRACKET m_material RIGHT_BRACKET SEMICOLON {
dynare_parser->set_vcov_pos(@4.off, @5.off);
};
m_material : m_material CHARACTER | m_material NAME | m_material SEMICOLON | m_material COMMA | CHARACTER | NAME | SEMICOLON | COMMA;
order : ORDER EQUAL_SIGN material SEMICOLON {
dynare_parser->set_order_pos(@3.off, @4.off);
};
planner : planner_objective planner_discount
| planner_discount planner_objective
;
planner_objective : PLANNEROBJECTIVE material SEMICOLON {
dynare_parser->set_pl_objective_pos(@2.off, @3.off);
};
planner_discount : PLANNERDISCOUNT NAME SEMICOLON {
dynare_parser->set_pl_discount_pos(@2.off, @3.off);
};
%%
void dynglob_error(const char* mes)
{
dynare_parser->error(mes);
}
/*
static void print_token_value1(FILE* file, int type, YYSTYPE value)
{
if (type == NAME)
fprintf(file, "%s", value.string);
if (type == CHARACTER)
fprintf(file, "%c", value.character);
}
*/
dynare++/src/forw_subst_builder.cpp deleted 100644 → 0
View file @ 0691b083
// Copyright (C) 2006-2011, Ondra Kamenik
#include "forw_subst_builder.h"
using namespace ogdyn;
ForwSubstBuilder::ForwSubstBuilder(DynareModel& m)
: model(m)
{
info.num_new_terms -= model.getParser().getTree().get_num_op();
// go through all equations
int neq = model.eqs.nformulas();
for (int i = 0; i < neq; i++) {
int ft = model.eqs.formula(i);
int mlead, mlag;
model.termspan(ft, mlead, mlag);
// if equation is too forward looking
if (mlead > 1) {
info.num_affected_equations++;
// break it to non-linear terms
unordered_set<int> nlt = model.get_nonlinear_subterms(ft);
int j = 0; // indexes subterms
// and make substitutions for all these non-linear subterms
for (unordered_set<int>::const_iterator it = nlt.begin();
it != nlt.end(); ++it, ++j)
substitute_for_term(*it, i, j);
}
}
// unassign all variables with lead greater than 1
unassign_gt_1_leads();
// forget the derivatives in the tree because some variables could
// have been unassigned
model.eqs.getTree().forget_derivative_maps();
info.num_new_terms += model.getParser().getTree().get_num_op();
}
void ForwSubstBuilder::substitute_for_term(int t, int i, int j)
{
int mlead, mlag;
model.termspan(t, mlead, mlag);
if (mlead > 1) {
info.num_subst_terms++;
// Example for comments: let t = f(x(+4))
// first make lagsubst be substitution setting f(x(+4)) to f(x(+1))
// this is lag = -3 (1-mlead)
map<int,int> lagsubst;
model.variable_shift_map(model.eqs.nulary_of_term(t), 1-mlead, lagsubst);
int lagt = model.eqs.add_substitution(t, lagsubst);
// now maxlead of lagt is +1
// add AUXLD_*_*_1 = f(x(+1)) to the model
char name[100];
sprintf(name, "AUXLD_%d_%d_%d", i, j, 1);
model.atoms.register_uniq_endo(name);
info.num_aux_variables++;
const char* ss = model.atoms.get_name_storage().query(name);
int auxt = model.eqs.add_nulary(name);
model.eqs.add_formula(model.eqs.add_binary(ogp::MINUS, auxt, lagt));
aux_map.insert(Tsubstmap::value_type(ss, lagt));
// now add variables and equations
// AUXLD_*_*_2 = AUXLD_*_*_1(+1) through
// AUXLD_*_*_{mlead-1} = AUXLD_*_*_{mlead-2}(+1)
for (int ll = 1; ll <= mlead-2; ll++) {
// create AUXLD_*_*_{ll}(+1)
sprintf(name, "AUXLD_%d_%d_%d(+1)", i, j, ll);
int lastauxt_lead = model.eqs.add_nulary(name);
// create AUXLD_*_*{ll+1}
sprintf(name, "AUXLD_%d_%d_%d", i, j, ll+1);
model.atoms.register_uniq_endo(name);
info.num_aux_variables++;
ss = model.atoms.get_name_storage().query(name);
auxt = model.eqs.add_nulary(name);
// add AUXLD_*_*_{ll+1} = AUXLD_*_*_{ll}(+1)
model.eqs.add_formula(model.eqs.add_binary(ogp::MINUS, auxt, lastauxt_lead));
// add substitution to the map; todo: this
// works well because in the context where
// aux_map is used the timing doesn't matter,
// however, it is misleading, needs to be
// changed
aux_map.insert(Tsubstmap::value_type(ss, lagt));
}
// now we have to substitute AUXLEAD_*_*{mlead-1}(+1) for t
model.substitute_atom_for_term(ss, +1, t);
}
}
void ForwSubstBuilder::unassign_gt_1_leads(const char* name)
{
const char* ss = model.atoms.get_name_storage().query(name);
int mlead, mlag;
model.atoms.varspan(name, mlead, mlag);
for (int ll = 2; ll <= mlead; ll++) {
int t = model.atoms.index(ss, ll);
if (t != -1)
model.atoms.unassign_variable(ss, ll, t);
}
}
void ForwSubstBuilder::unassign_gt_1_leads()
{
const vector<const char*>& endovars = model.atoms.get_endovars();
for (unsigned int i = 0; i < endovars.size(); i++)
unassign_gt_1_leads(endovars[i]);
const vector<const char*>& exovars = model.atoms.get_exovars();
for (unsigned int i = 0; i < exovars.size(); i++)
unassign_gt_1_leads(exovars[i]);
}
ForwSubstBuilder::ForwSubstBuilder(const ForwSubstBuilder& b, DynareModel& m)
: model(m)
{
for (Tsubstmap::const_iterator it = b.aux_map.begin();
it != b.aux_map.end(); ++it) {
const char* ss = m.atoms.get_name_storage().query((*it).first);
aux_map.insert(Tsubstmap::value_type(ss, (*it).second));
}
}
dynare++/src/forw_subst_builder.h deleted 100644 → 0
View file @ 0691b083
// Copyright (C) 2006, Ondra Kamenik
// $Id$
#ifndef FORW_SUBST_BUILDER_H
#define FORW_SUBST_BUILDER_H
#include "dynare_model.h"
namespace ogdyn {
/** This struct encapsulates information about the process of
* forward substitutions. */
struct ForwSubstInfo {
int num_affected_equations;
int num_subst_terms;
int num_aux_variables;
int num_new_terms;
ForwSubstInfo()
: num_affected_equations(0),
num_subst_terms(0),
num_aux_variables(0),
num_new_terms(0) {}
};
class ForwSubstBuilder {
typedef map<int, const char*> Ttermauxmap;
protected:
/** Reference to the model, to which we will add equations and
* change some equations. */
DynareModel& model;
/** A map mapping new auxiliary variables to the terms in the
* tree in the DynareModel. */
Tsubstmap aux_map;
/** Information about the substitutions. */
ForwSubstInfo info;
public:
/** Do all the jobs needed. This scans all equations in the
* model, and for equations containing forward looking
* variables greater than 1 lead, it makes corresponding
* substitutions. Basically, it breaks each equation to its
* non-linear components and creates substitutions for these
* components, not for whole equation. This is because the
* expectation operator can go through the linear part of the
* function. This will save us many occurrences of other
* variables involved in the equation. */
ForwSubstBuilder(DynareModel& m);
/** Copy constructor with a new instance of the model. */
ForwSubstBuilder(const ForwSubstBuilder& b, DynareModel& m);
/** Return the auxiliary variable mapping. */
const Tsubstmap& get_aux_map() const
{return aux_map;}
/** Return the information. */
const ForwSubstInfo& get_info() const
{return info;}
private:
ForwSubstBuilder(const ForwSubstBuilder& b);
/** This method takes a nonlinear term t, and if it has leads
* of greater than 1, then it substitutes the term for the new
* variable (or string of variables). Note that the
* substitution is done by DynamicAtoms::assign_variable. This
* means that the substitution is made for all other
* ocurrences of t in the model. So there is no need of
* tracking already substituted terms. The other two
* parameters are just for identification of the new auxiliary
* variables. When called from the constructor, i is an
* equation number, j is an order of the non-linear term in
* the equation. */
void substitute_for_term(int t, int i, int j);
/** This is called just at the end of the job. It unassigns
* all nulary terms with a lead greater than 1. */
void unassign_gt_1_leads();
/** This unassigns all leads greater than 1 of the given name. */
void unassign_gt_1_leads(const char* name);
};
};
#endif
// Local Variables:
// mode:C++
// End:
dynare++/src/main.cpp deleted 100644 → 0
View file @ 0691b083
// Copyright (C) 2004-2011, Ondra Kamenik
#include "dynare3.h"
#include "dynare_exception.h"
#include "dynare_params.h"
#include "utils/cc/exception.h"
#include "parser/cc/parser_exception.h"
#include "../sylv/cc/SylvException.h"
#include "../kord/random.h"
#include "../kord/global_check.h"
#include "../kord/approximation.h"
int main(int argc, char** argv)
{
DynareParams params(argc, argv);
if (params.help) {
params.printHelp();
return 0;
}
if (params.version) {
printf("Dynare++ v. %s. Copyright (C) 2004-2011, Ondra Kamenik\n",
DYNVERSION);
printf("Dynare++ comes with ABSOLUTELY NO WARRANTY and is distributed under\n");
printf("GPL: modules integ, tl, kord, sylv, src, extern and documentation\n");
printf("LGPL: modules parser, utils\n");
printf(" for GPL see http://www.gnu.org/licenses/gpl.html\n");
printf(" for LGPL see http://www.gnu.org/licenses/lgpl.html\n");
return 0;
}
THREAD_GROUP::max_parallel_threads = params.num_threads;
try {
// make journal name and journal
std::string jname(params.basename);
jname += ".jnl";
Journal journal(jname.c_str());
// make dynare object
Dynare dynare(params.modname, params.order, params.ss_tol, journal);
// make list of shocks for which we will do IRFs
vector<int> irf_list_ind;
if (params.do_irfs_all)
for (int i = 0; i < dynare.nexog(); i++)
irf_list_ind.push_back(i);
else
irf_list_ind = ((const DynareNameList&)dynare.getExogNames()).selectIndices(params.irf_list);
// write matlab files
FILE* mfd;
std::string mfile1(params.basename);
mfile1 += "_f.m";
if (NULL == (mfd=fopen(mfile1.c_str(), "w"))) {
fprintf(stderr, "Couldn't open %s for writing.\n", mfile1.c_str());
exit(1);
}
ogdyn::MatlabSSWriter writer0(dynare.getModel(), params.basename.c_str());
writer0.write_der0(mfd);
fclose(mfd);
std::string mfile2(params.basename);
mfile2 += "_ff.m";
if (NULL == (mfd=fopen(mfile2.c_str(), "w"))) {
fprintf(stderr, "Couldn't open %s for writing.\n", mfile2.c_str());
exit(1);
}
ogdyn::MatlabSSWriter writer1(dynare.getModel(), params.basename.c_str());
writer1.write_der1(mfd);
fclose(mfd);
// open mat file
std::string matfile(params.basename);
matfile += ".mat";
mat_t* matfd = Mat_Create(matfile.c_str(), NULL);
if (matfd == NULL) {
fprintf(stderr, "Couldn't open %s for writing.\n", matfile.c_str());
exit(1);
}
// write info about the model (dimensions and variables)
dynare.writeMat(matfd, params.prefix);
// write the dump file corresponding to the input
dynare.writeDump(params.basename);
system_random_generator.initSeed(params.seed);
tls.init(dynare.order(),
dynare.nstat()+2*dynare.npred()+3*dynare.nboth()+
2*dynare.nforw()+dynare.nexog());
Approximation app(dynare, journal, params.num_steps, params.do_centralize, params.qz_criterium);
try {
app.walkStochSteady();
} catch (const KordException& e) {
// tell about the exception and continue
printf("Caught (not yet fatal) Kord exception: ");
e.print();
JournalRecord rec(journal);
rec << "Solution routine not finished (" << e.get_message()
<< "), see what happens" << endrec;
}
std::string ss_matrix_name(params.prefix);
ss_matrix_name += "_steady_states";
ConstTwoDMatrix(app.getSS()).writeMat(matfd, ss_matrix_name.c_str());
// check the approximation
if (params.check_along_path || params.check_along_shocks
|| params.check_on_ellipse) {
GlobalChecker gcheck(app, THREAD_GROUP::max_parallel_threads, journal);
if (params.check_along_shocks)
gcheck.checkAlongShocksAndSave(matfd, params.prefix,
params.getCheckShockPoints(),
params.getCheckShockScale(),
params.check_evals);
if (params.check_on_ellipse)
gcheck.checkOnEllipseAndSave(matfd, params.prefix,
params.getCheckEllipsePoints(),
params.getCheckEllipseScale(),
params.check_evals);
if (params.check_along_path)
gcheck.checkAlongSimulationAndSave(matfd, params.prefix,
params.getCheckPathPoints(),
params.check_evals);
}
// write the folded decision rule to the Mat-4 file
app.getFoldDecisionRule().writeMat(matfd, params.prefix);
// simulate conditional
if (params.num_condper > 0 && params.num_condsim > 0) {
SimResultsDynamicStats rescond(dynare.numeq(), params.num_condper, 0);
ConstVector det_ss(app.getSS(),0);
rescond.simulate(params.num_condsim, app.getFoldDecisionRule(), det_ss, dynare.getVcov(), journal);
rescond.writeMat(matfd, params.prefix);
}
// simulate unconditional
//const DecisionRule& dr = app.getUnfoldDecisionRule();
const DecisionRule& dr = app.getFoldDecisionRule();
if (params.num_per > 0 && params.num_sim > 0) {
SimResultsStats res(dynare.numeq(), params.num_per, params.num_burn);
res.simulate(params.num_sim, dr, dynare.getSteady(), dynare.getVcov(), journal);
res.writeMat(matfd, params.prefix);
// impulse response functions
if (! irf_list_ind.empty()) {
IRFResults irf(dynare, dr, res, irf_list_ind, journal);
irf.writeMat(matfd, params.prefix);
}
}
// simulate with real-time statistics
if (params.num_rtper > 0 && params.num_rtsim > 0) {
RTSimResultsStats rtres(dynare.numeq(), params.num_rtper, params.num_burn);
rtres.simulate(params.num_rtsim, dr, dynare.getSteady(), dynare.getVcov(), journal);
rtres.writeMat(matfd, params.prefix);
}
Mat_Close(matfd);
} catch (const KordException& e) {
printf("Caugth Kord exception: ");
e.print();
return e.code();
} catch (const TLException& e) {
printf("Caugth TL exception: ");
e.print();
return 255;
} catch (SylvException& e) {
printf("Caught Sylv exception: ");
e.printMessage();
return 255;
} catch (const DynareException& e) {
printf("Caught Dynare exception: %s\n", e.message());
return 255;
} catch (const ogu::Exception& e) {
printf("Caught ogu::Exception: ");
e.print();
return 255;
} catch (const ogp::ParserException& e) {
printf("Caught parser exception: %s\n", e.message());
return 255;
}
return 0;
}
dynare++/src/nlsolve.cpp deleted 100644 → 0
View file @ 0691b083
// Copyright (C) 2006, Ondra Kamenik
// $Id: nlsolve.cpp 762 2006-05-22 13:00:07Z kamenik $
#include "nlsolve.h"
#include "dynare_exception.h"
#include <cmath>
using namespace ogu;
/** This should not be greater than DBL_EPSILON^(1/2). */
double GoldenSectionSearch::tol = 1.e-4;
/** This is equal to the golden section ratio. */
double GoldenSectionSearch::golden = (3.-std::sqrt(5.))/2;
double GoldenSectionSearch::search(OneDFunction& f, double x1, double x2)
{
double b;
if (init_bracket(f, x1, x2, b)) {
double fb = f.eval(b);
double f1 = f.eval(x1);
f.eval(x2);
double dx;
do {
double w = (b-x1)/(x2-x1);
dx = std::abs((1-2*w)*(x2-x1));
double x;
if (b-x1 > x2-b)
x = b - dx;
else
x = b + dx;
double fx = f.eval(x);
if (! std::isfinite(fx))
return x1;
if (b-x1 > x2-b) {
// x is on the left from b
if (f1 > fx && fx < fb) {
// pickup bracket [f1,fx,fb]
x2 = b;
fb = fx;
b = x;
} else {
// pickup bracket [fx,fb,fx2]
f1 = fx;
x1 = x;
}
} else {
// x is on the right from b
if (f1 > fb && fb < fx) {
// pickup bracket [f1,fb,fx]
x2 = x;
} else {
// pickup bracket [fb,fx,f2]
f1 = fb;
x1 = b;
fb = fx;
b = x;
}
}
} while(dx > tol);
}
return b;
}
bool GoldenSectionSearch::init_bracket(OneDFunction& f, double x1, double& x2, double& b)
{
double f1 = f.eval(x1);
if (! std::isfinite(f1))
throw DynareException(__FILE__, __LINE__,
"Safer point not finite in GoldenSectionSearch::init_bracket");
int cnt = 0;
bool bracket_found = false;
do {
bool finite_found = search_for_finite(f, x1, x2, b);
if (! finite_found) {
b = x1;
return false;
}
double f2 = f.eval(x2);
double fb = f.eval(b);
double bsym = 2*x2 - b;
double fbsym = f.eval(bsym);
// now we know that f1, f2, and fb are finite
if (std::isfinite(fbsym)) {
// we have four numbers f1, fb, f2, fbsym, we test for the
// following combinations to find the bracket:
// [f1,f2,fbsym], [f1,fb,fbsym] and [f1,fb,fbsym]
if (f1 > f2 && f2 < fbsym) {
bracket_found = true;
b = x2;
x2 = bsym;
} else if (f1 > fb && fb < fbsym) {
bracket_found = true;
x2 = bsym;
} else if (f1 > fb && fb < f2) {
bracket_found = true;
} else {
double newx2 = b;
// choose the smallest value in case we end
if (f1 > fbsym) {
// the smallest value is on the other end, we do
// not want to continue
b = bsym;
return false;
} else
b = x1;
// move x2 to b in case we continue
x2 = newx2;
}
} else {
// we have only three numbers, we test for the bracket,
// and if not found, we set b as potential result and
// shorten x2 as potential init value for next cycle
if (f1 > fb && fb < f2)
bracket_found = true;
else {
double newx2 = b;
// choose the smaller value in case we end
if (f1 > f2)
b = x2;
else
b = x1;
// move x2 to b in case we continue
x2 = newx2;
}
}
cnt++;
} while (! bracket_found && cnt < 5);
return bracket_found;
}
/** This moves x2 toward to x1 until the function at x2 is finite and
* b as a golden section between x1 and x2 yields also finite f. */
bool GoldenSectionSearch::search_for_finite(OneDFunction& f, double x1, double& x2, double&b)
{
int cnt = 0;
bool found = false;
do {
double f2 = f.eval(x2);
b = (1-golden)*x1 + golden*x2;
double fb = f.eval(b);
found = std::isfinite(f2) && std::isfinite(fb);
if (! found)
x2 = b;
cnt++;
} while (! found && cnt < 5);
return found;
}
void VectorFunction::check_for_eval(const ConstVector& in, Vector& out) const
{
if (inDim() != in.length() || outDim() != out.length())
throw DynareException(__FILE__, __LINE__,
"Wrong dimensions in VectorFunction::check_for_eval");
}
double NLSolver::eval(double lambda)
{
Vector xx((const Vector&)x);
xx.add(1-lambda, xcauchy);
xx.add(lambda, xnewton);
Vector ff(func.outDim());
func.eval(xx, ff);
return ff.dot(ff);
}
bool NLSolver::solve(Vector& xx, int& iter)
{
JournalRecord rec(journal);
rec << "Iter lambda residual" << endrec;
JournalRecord rec1(journal);
rec1 << "---------------------------" << endrec;
char tmpbuf[14];
x = (const Vector&)xx;
iter = 0;
// setup fx
Vector fx(func.outDim());
func.eval(x, fx);
if (!fx.isFinite())
throw DynareException(__FILE__,__LINE__,
"Initial guess does not yield finite residual in NLSolver::solve");
bool converged = fx.getMax() < tol;
JournalRecord rec2(journal);
sprintf(tmpbuf, "%10.6g", fx.getMax());
rec2 << iter << " N/A " << tmpbuf << endrec;
while (! converged && iter < max_iter) {
// setup Jacobian
jacob.eval(x);
// calculate cauchy step
Vector g(func.inDim());
g.zeros();
ConstTwoDMatrix(jacob).multaVecTrans(g, fx);
Vector Jg(func.inDim());
Jg.zeros();
ConstTwoDMatrix(jacob).multaVec(Jg, g);
double m = -g.dot(g)/Jg.dot(Jg);
xcauchy = (const Vector&) g;
xcauchy.mult(m);
// calculate newton step
xnewton = (const Vector&) fx;
ConstTwoDMatrix(jacob).multInvLeft(xnewton);
xnewton.mult(-1);
// line search
double lambda = GoldenSectionSearch::search(*this, 0, 1);
x.add(1-lambda, xcauchy);
x.add(lambda, xnewton);
// evaluate func
func.eval(x, fx);
converged = fx.getMax() < tol;
// iter
iter++;
JournalRecord rec3(journal);
sprintf(tmpbuf, "%10.6g", fx.getMax());
rec3 << iter << " " << lambda << " " << tmpbuf << endrec;
}
xx = (const Vector&)x;
return converged;
}
dynare++/src/nlsolve.h deleted 100644 → 0
View file @ 0691b083
// Copyright (C) 2006, Ondra Kamenik
// $Id: nlsolve.h 762 2006-05-22 13:00:07Z kamenik $
#ifndef OGU_NLSOLVE_H
#define OGU_NLSOLVE_H
#include "twod_matrix.h"
#include "journal.h"
namespace ogu {
class OneDFunction {
public:
virtual ~OneDFunction() {}
virtual double eval(double) = 0;
};
class GoldenSectionSearch {
protected:
static double tol;
static double golden;
public:
static double search(OneDFunction& f, double x1, double x2);
protected:
/** This initializes a bracket by moving x2 and b (as a golden
* section of x1,x2) so that f(x1)>f(b) && f(b)<f(x2). The point
* x1 is not moved, since it is considered as reliable and f(x1)
* is supposed to be finite. If initialization of a bracket
* succeeded, then [x1, b, x2] is the bracket and true is
* returned. Otherwise, b is the minimum found and false is
* returned. */
static bool init_bracket(OneDFunction& f, double x1, double& x2, double& b);
/** This supposes that f(x1) is finite and it moves x2 toward x1
* until x2 and b (as a golden section of x1,x2) are finite. If
* succeeded, the routine returns true and x2, and b. Otherwise,
* it returns false. */
static bool search_for_finite(OneDFunction& f, double x1, double& x2, double& b);
};
class VectorFunction {
public:
VectorFunction() {}
virtual ~VectorFunction() {}
virtual int inDim() const = 0;
virtual int outDim() const = 0;
/** Check dimensions of eval parameters. */
void check_for_eval(const ConstVector& in, Vector& out) const;
/** Evaluate the vector function. */
virtual void eval(const ConstVector& in, Vector& out) = 0;
};
class Jacobian : public TwoDMatrix {
public:
Jacobian(int n)
: TwoDMatrix(n,n) {}
virtual ~Jacobian() {}
virtual void eval(const Vector& in) = 0;
};
class NLSolver : public OneDFunction {
protected:
Journal& journal;
VectorFunction& func;
Jacobian& jacob;
const int max_iter;
const double tol;
private:
Vector xnewton;
Vector xcauchy;
Vector x;
public:
NLSolver(VectorFunction& f, Jacobian& j, int maxit, double tl, Journal& jr)
: journal(jr), func(f), jacob(j), max_iter(maxit), tol(tl),
xnewton(f.inDim()), xcauchy(f.inDim()), x(f.inDim())
{xnewton.zeros(); xcauchy.zeros(); x.zeros();}
virtual ~NLSolver() {}
/** Returns true if the problem has converged. xx as input is the
* starting value, as output it is a solution. */
bool solve(Vector& xx, int& iter);
/** To implement OneDFunction interface. It returns
* func(xx)^T*func(xx), where
* xx=x+lambda*xcauchy+(1-lambda)*xnewton. It is non-const only
* because it calls func, x, xnewton, xcauchy is not changed. */
double eval(double lambda);
};
};
#endif
// Local Variables:
// mode:C++
// End:
dynare++/src/planner_builder.cpp deleted 100644 → 0
View file @ 0691b083
// Copyright (C) 2006, Ondra Kamenik
// $Id$
#include "planner_builder.h"
#include "dynare_exception.h"
#include <cmath>
using namespace ogdyn;
const IntegerMatrix& IntegerMatrix::operator=(const IntegerMatrix& im)
{
if (nr != im.nr || nc != im.nc)
throw DynareException(__FILE__,__LINE__,
"Matrices have different dimensions in IntegerMatrix::operator=");
memcpy(data, im.data, nr*nc*sizeof(int));
return *this;
}
const IntegerArray3& IntegerArray3::operator=(const IntegerArray3& ia3)
{
if (n1 != ia3.n1 || n2 != ia3.n2 || n3 != ia3.n3)
throw DynareException(__FILE__,__LINE__,
"Arrays have different dimensions in IntegerArray3::operator=");
memcpy(data, ia3.data, n1*n2*n3*sizeof(int));
return *this;
}
PlannerBuilder::PlannerBuilder(DynareModel& m, const Tvarset& yyset,
const Teqset& ffset)
: yset(), fset(ffset), model(m),
tb(model.t_plobjective), tbeta(model.t_pldiscount),
maxlead(model.atoms.get_maxlead()),
minlag(model.atoms.get_minlag()),
diff_b(yyset.size(), 1-minlag),
diff_f(yyset.size(), fset.size(), 1+maxlead-minlag),
static_atoms(),
static_tree(),
diff_b_static(yyset.size(), 1-minlag),
diff_f_static(yyset.size(), fset.size(), 1+maxlead-minlag)
{
info.num_new_terms -= model.getParser().getTree().get_num_op();
fill_yset(m.atoms.get_name_storage(), yyset);
add_derivatives_of_b();
add_derivatives_of_f();
shift_derivatives_of_b();
shift_derivatives_of_f();
beta_multiply_b();
beta_multiply_f();
make_static_version();
lagrange_mult_f();
form_equations();
info.num_new_terms += model.getParser().getTree().get_num_op();
}
PlannerBuilder::PlannerBuilder(const PlannerBuilder& pb, ogdyn::DynareModel& m)
: yset(), fset(pb.fset), model(m),
tb(pb.tb), tbeta(pb.tbeta),
maxlead(pb.maxlead), minlag(pb.minlag),
diff_b(pb.diff_b), diff_f(pb.diff_f),
static_atoms(pb.static_atoms),
static_tree(pb.static_tree),
diff_b_static(pb.diff_b_static),
diff_f_static(pb.diff_f_static),
aux_map(), static_aux_map()
{
fill_yset(m.atoms.get_name_storage(), pb.yset);
fill_aux_map(m.atoms.get_name_storage(), pb.aux_map, pb.static_aux_map);
}
void PlannerBuilder::add_derivatives_of_b()
{
int yi = 0;
for (Tvarset::const_iterator yname = yset.begin();
yname != yset.end(); ++yname, yi++)
for (int ll = minlag; ll <= 0; ll++) {
int yt = model.atoms.index(*yname, ll);
if (yt != -1)
diff_b(yi, ll-minlag) = model.eqs.add_derivative(tb, yt);
else
diff_b(yi, ll-minlag) = ogp::OperationTree::zero;
}
}
void PlannerBuilder::add_derivatives_of_f()
{
int yi = 0;
for (Tvarset::const_iterator yname = yset.begin();
yname != yset.end(); ++yname, yi++)
for (unsigned int fi = 0; fi < fset.size(); fi++)
for (int ll = minlag; ll <= maxlead; ll++) {
int yt = model.atoms.index(*yname, ll);
if (yt != -1)
diff_f(yi, fi, ll-minlag) =
model.eqs.add_derivative(model.eqs.formula(fset[fi]), yt);
else
diff_f(yi, fi, ll-minlag) = ogp::OperationTree::zero;
}
}
void PlannerBuilder::shift_derivatives_of_b()
{
map<int,int> subst;
for (int yi = 0; yi < diff_b.nrows(); yi++)
for (int ll = minlag; ll < 0; ll++)
if (diff_b(yi, ll-minlag) != ogp::OperationTree::zero) {
model.variable_shift_map(model.eqs.nulary_of_term(diff_b(yi, ll-minlag)),
-ll, subst);
diff_b(yi, ll-minlag) = model.eqs.add_substitution(diff_b(yi, ll-minlag), subst);
}
}
void PlannerBuilder::shift_derivatives_of_f()
{
map<int,int> subst;
for (int yi = 0; yi < diff_f.dim1(); yi++)
for (int fi = 0; fi < diff_f.dim2(); fi++) {
// first do it leads which are put under expectation before t: no problem
for (int ll = 0; ll <= maxlead; ll++)
if (diff_f(yi, fi, ll-minlag) != ogp::OperationTree::zero) {
model.variable_shift_map(model.eqs.nulary_of_term(diff_f(yi, fi, ll-minlag)),
-ll, subst);
diff_f(yi, fi, ll-minlag) =
model.eqs.add_substitution(diff_f(yi, fi, ll-minlag), subst);
}
// now do it for lags, these are put as leads under
// expectations after time t, so we have to introduce
// auxiliary variables at time t, and make leads of them here
for (int ll = minlag; ll < 0; ll++) {
int ft = diff_f(yi, fi, ll-minlag);
if (ft != ogp::OperationTree::zero) {
// if the ft term has a lead, than we need to
// introduce an auxiliary variable z_t, define it
// as E_t[ft] and put z_{t-ll} to the
// equation. Otherwise, we just put leaded ft to
// the equation directly.
int ft_maxlead, ft_minlag;
model.termspan(ft, ft_maxlead, ft_minlag);
if (ft_maxlead > 0) {
// make an auxiliary variable
char name[100];
sprintf(name, "AUX_%d_%d_%d", yi, fset[fi], -ll);
model.atoms.register_uniq_endo(name);
info.num_aux_variables++;
int taux = model.eqs.add_nulary(name);
sprintf(name, "AUX_%d_%d_%d(%d)", yi, fset[fi], -ll, -ll);
int taux_leaded = model.eqs.add_nulary(name);
// put aux_leaded to the equation
diff_f(yi, fi, ll-minlag) = taux_leaded;
// save auxiliary variable and the term
aux_map.insert(Tsubstmap::value_type(model.atoms.name(taux), ft));
} else {
// no auxiliary variable is needed and the
// term ft can be leaded in place
model.variable_shift_map(model.eqs.nulary_of_term(ft), -ll, subst);
diff_f(yi, fi, ll-minlag) =
model.eqs.add_substitution(ft, subst);
}
}
}
}
}
void PlannerBuilder::beta_multiply_b()
{
int beta_pow = ogp::OperationTree::one;
for (int ll = 0; ll >= minlag; ll--,
beta_pow = model.eqs.add_binary(ogp::TIMES, beta_pow, tbeta))
for (int yi = 0; yi < diff_b.nrows(); yi++)
if (diff_b(yi, ll-minlag) != ogp::OperationTree::zero)
diff_b(yi, ll-minlag) =
model.eqs.add_binary(ogp::TIMES, beta_pow, diff_b(yi, ll-minlag));
}
void PlannerBuilder::beta_multiply_f()
{
int beta_pow = ogp::OperationTree::one;
for (int ll = 0; ll <= maxlead; ll++,
beta_pow = model.eqs.add_binary(ogp::DIVIDE, beta_pow, tbeta))
for (int yi = 0; yi < diff_f.dim1(); yi++)
for (int fi = 0; fi < diff_f.dim2(); fi++)
if (diff_f(yi, fi, ll-minlag) != ogp::OperationTree::zero)
diff_f(yi, fi, ll-minlag) =
model.eqs.add_binary(ogp::TIMES, beta_pow, diff_f(yi, fi, ll-minlag));
beta_pow = ogp::OperationTree::one;
for (int ll = 0; ll >= minlag; ll--,
beta_pow = model.eqs.add_binary(ogp::TIMES, beta_pow, tbeta))
for (int yi = 0; yi < diff_f.dim1(); yi++)
for (int fi = 0; fi < diff_f.dim2(); fi++)
if (diff_f(yi, fi, ll-minlag) != ogp::OperationTree::zero)
diff_f(yi, fi, ll-minlag) =
model.eqs.add_binary(ogp::TIMES, beta_pow, diff_f(yi, fi, ll-minlag));
}
void PlannerBuilder::make_static_version()
{
// map holding substitutions from dynamic to static
ogp::StaticFineAtoms::Tintintmap tmap;
// fill static atoms with outer ordering
static_atoms.import_atoms(model.atoms, static_tree, tmap);
// go through diff_b and fill diff_b_static
for (int ll = minlag; ll <= 0; ll++)
for (int yi = 0; yi < diff_b.nrows(); yi++)
diff_b_static(yi, ll-minlag) =
static_tree.add_substitution(diff_b(yi, ll-minlag),
tmap, model.eqs.getTree());
// go through diff_f and fill diff_f_static
for (int ll = minlag; ll <= maxlead; ll++)
for (int yi = 0; yi < diff_f.dim1(); yi++)
for (int fi = 0; fi < diff_f.dim2(); fi++)
diff_f_static(yi, fi, ll-minlag) =
static_tree.add_substitution(diff_f(yi, fi, ll-minlag),
tmap, model.eqs.getTree());
// go through aux_map and fill static_aux_map
for (Tsubstmap::const_iterator it = aux_map.begin();
it != aux_map.end(); ++it) {
int tstatic = static_tree.add_substitution((*it).second, tmap, model.eqs.getTree());
const char* name = static_atoms.get_name_storage().query((*it).first);
static_aux_map.insert(Tsubstmap::value_type(name, tstatic));
}
}
void PlannerBuilder::lagrange_mult_f()
{
// register multipliers
char mult_name[100];
for (int fi = 0; fi < diff_f.dim2(); fi++) {
sprintf(mult_name, "MULT%d", fset[fi]);
model.atoms.register_uniq_endo(mult_name);
info.num_lagrange_mults++;
}
// multiply with the multipliers
for (int yi = 0; yi < diff_f.dim1(); yi++)
for (int fi = 0; fi < diff_f.dim2(); fi++)
for (int ll = minlag; ll <= maxlead; ll++)
if (diff_f(yi, fi, ll-minlag) != ogp::OperationTree::zero) {
sprintf(mult_name, "MULT%d(%d)", fset[fi], -ll);
int tm = model.eqs.add_nulary(mult_name);
diff_f(yi, fi, ll-minlag) =
model.eqs.add_binary(ogp::TIMES, tm, diff_f(yi, fi, ll-minlag));
}
}
void PlannerBuilder::form_equations()
{
// add planner's FOCs
for (int yi = 0; yi < diff_f.dim1(); yi++) {
int eq = ogp::OperationTree::zero;
for (int ll = minlag; ll <= 0; ll++)
eq = model.eqs.add_binary(ogp::PLUS, eq, diff_b(yi, ll-minlag));
for (int fi = 0; fi < diff_f.dim2(); fi++)
for (int ll = minlag; ll <= maxlead; ll++)
eq = model.eqs.add_binary(ogp::PLUS, eq, diff_f(yi, fi, ll-minlag));
model.eqs.add_formula(eq);
}
// add equations for auxiliary variables
for (Tsubstmap::const_iterator it = aux_map.begin();
it != aux_map.end(); ++it) {
int t = model.atoms.index((*it).first, 0);
model.eqs.add_formula(model.eqs.add_binary(ogp::MINUS, t, (*it).second));
}
}
void PlannerBuilder::fill_yset(const ogp::NameStorage& ns,
const PlannerBuilder::Tvarset& yyset)
{
for (Tvarset::const_iterator it = yyset.begin(); it != yyset.end(); ++it)
yset.insert(ns.query(*it));
}
void PlannerBuilder::fill_aux_map(const ogp::NameStorage& ns, const Tsubstmap& aaux_map,
const Tsubstmap& astatic_aux_map)
{
// fill aux_map
for (Tsubstmap::const_iterator it = aaux_map.begin();
it != aaux_map.end(); ++it)
aux_map.insert(Tsubstmap::value_type(ns.query((*it).first), (*it).second));
// fill static_aux_map
for (Tsubstmap::const_iterator it = astatic_aux_map.begin();
it != astatic_aux_map.end(); ++it)
static_aux_map.insert(Tsubstmap::value_type(static_atoms.get_name_storage().query((*it).first),
(*it).second));
}
MultInitSS::MultInitSS(const PlannerBuilder& pb, const Vector& pvals, Vector& yy)
: builder(pb), b(builder.diff_b_static.nrows()),
F(builder.diff_f_static.dim1(), builder.diff_f_static.dim2())
{
b.zeros();
F.zeros();
// first evaluate substitutions (auxiliary variables) from the builder
ogdyn::DynareStaticSteadySubstitutions dss(builder.model.atoms, builder.static_atoms,
builder.static_tree,
builder.static_aux_map, pvals, yy);
// gather all the terms from builder.diff_b_static and
// builder.diff_f_static to the vector, the ordering is important,
// since the index of this vector will have to be decoded to the
// position in b and F.
vector<int> terms;
for (int yi = 0; yi < builder.diff_b_static.nrows(); yi++)
for (int l = 0; l < builder.diff_b_static.ncols(); l++)
terms.push_back(builder.diff_b_static(yi, l));
for (int yi = 0; yi < builder.diff_f_static.dim1(); yi++)
for (int fi = 0; fi < builder.diff_f_static.dim2(); fi++)
for (int l = 0; l < builder.diff_f_static.dim3(); l++)
terms.push_back(builder.diff_f_static(yi, fi, l));
// evaluate the terms, it will call a series of load(i,res), which
// sum the results through lags/leads to b and F
DynareStaticSteadyAtomValues dssav(builder.model.atoms, builder.static_atoms, pvals, yy);
ogp::FormulaCustomEvaluator fe(builder.static_tree, terms);
fe.eval(dssav, *this);
// solve overdetermined system b+F*lambda=0 => lambda=-(F^T*F)^{-1}*F^T*b
GeneralMatrix FtF(F, "transpose", F);
Vector lambda(builder.diff_f_static.dim2());
F.multVecTrans(0.0, lambda, -1.0, b);
ConstGeneralMatrix(FtF).multInvLeft(lambda);
// take values of lambda and put it to yy
for (int fi = 0; fi < builder.diff_f_static.dim2(); fi++) {
char mult_name[100];
sprintf(mult_name, "MULT%d", builder.fset[fi]);
int iouter = builder.model.atoms.name2outer_endo(mult_name);
int iy = builder.model.atoms.outer2y_endo()[iouter];
if (! std::isfinite(yy[iy]))
yy[iy] = lambda[fi];
// go through all substitutions of the multiplier and set them
// as well
if (builder.model.atom_substs) {
const ogp::AtomSubstitutions::Toldnamemap& old2new =
builder.model.atom_substs->get_old2new();
const ogp::AtomSubstitutions::Toldnamemap::const_iterator it =
old2new.find(mult_name);
if (it != old2new.end()) {
const ogp::AtomSubstitutions::Tshiftnameset& sset = (*it).second;
for (ogp::AtomSubstitutions::Tshiftnameset::const_iterator itt = sset.begin();
itt != sset.end(); ++itt) {
const char* newname = (*itt).first;
int iouter = builder.model.atoms.name2outer_endo(newname);
int iy = builder.model.atoms.outer2y_endo()[iouter];
if (! std::isfinite(yy[iy]))
yy[iy] = lambda[fi];
}
}
}
}
}
void MultInitSS::load(int i, double res)
{
// we can afford it, since the evaluator sets res to exact zero if
// the term is zero
if (res == 0)
return;
// decode i and add to either b or F
if (i < builder.diff_b_static.nrows()*builder.diff_b_static.ncols()) {
// add to b
b[i / builder.diff_b_static.ncols()] += res;
} else {
// add to F
i -= builder.diff_b_static.nrows()*builder.diff_b_static.ncols();
int yifi = i / builder.diff_f_static.dim3();
int yi = yifi / builder.diff_f_static.dim2();
int fi = yifi % builder.diff_f_static.dim2();
F.get(yi, fi) += res;
}
}
dynare++/src/planner_builder.h deleted 100644 → 0
View file @ 0691b083
// Copyright (C) 2006-2011, Ondra Kamenik
#ifndef PLANNER_BUILDER_H
#define PLANNER_BUILDER_H
#include "dynare_model.h"
namespace ogdyn {
using boost::unordered_set;
using std::map;
using std::vector;
/** This is a two dimensional array of integers. Nothing
* difficult. */
class IntegerMatrix {
protected:
/** Number of rows. */
int nr;
/** Number of columns. */
int nc;
/** The pointer to the data. */
int* data;
public:
/** Construct uninitialized array. */
IntegerMatrix(int nrr, int ncc)
: nr(nrr), nc(ncc), data(new int[nr*nc]) {}
/** Copy constructor. */
IntegerMatrix(const IntegerMatrix& im)
: nr(im.nr), nc(im.nc), data(new int[nr*nc])
{memcpy(data, im.data, nr*nc*sizeof(int));}
virtual ~IntegerMatrix()
{delete [] data;}
/** Assignment operator. It can only assing array with the
* same dimensions. */
const IntegerMatrix& operator=(const IntegerMatrix& im);
int& operator()(int i, int j)
{return data[i+j*nr];}
const int& operator()(int i, int j) const
{return data[i+j*nr];}
int nrows() const
{return nr;}
int ncols() const
{return nc;}
};
/** The three dimensional array of integers. Nothing difficult. */
class IntegerArray3 {
protected:
/** First dimension. */
int n1;
/** Second dimension. */
int n2;
/** Third dimension. */
int n3;
/** The data. */
int* data;
public:
/** Constrcut unitialized array. */
IntegerArray3(int nn1, int nn2, int nn3)
: n1(nn1), n2(nn2), n3(nn3), data(new int[n1*n2*n3]) {}
/** Copy constructor. */
IntegerArray3(const IntegerArray3& ia3)
: n1(ia3.n1), n2(ia3.n2), n3(ia3.n3), data(new int[n1*n2*n3])
{memcpy(data, ia3.data, n1*n2*n3*sizeof(int));}
virtual ~IntegerArray3()
{delete [] data;}
/** Assignment operator assigning the arrays with the same dimensions. */
const IntegerArray3& operator=(const IntegerArray3& ia3);
int& operator()(int i, int j, int k)
{return data[i+j*n1+k*n1*n2];}
const int& operator()(int i, int j, int k) const
{return data[i+j*n1+k*n1*n2];}
int dim1() const
{return n1;}
int dim2() const
{return n2;}
int dim3() const
{return n3;}
};
/** This struct encapsulates information about the building of a
* planner's problem. */
struct PlannerInfo {
int num_lagrange_mults;
int num_aux_variables;
int num_new_terms;
PlannerInfo()
: num_lagrange_mults(0),
num_aux_variables(0),
num_new_terms(0) {}
};
class MultInitSS;
/** This class builds the first order conditions of the social
* planner problem with constraints being the equations in the
* model. The model is non-const parameter to the constructor
* which adds appropriate FOCs to the system. It also allows for
* an estimation of the lagrange multipliers given all other
* endogenous variables of the static system. For this purpose we
* need to create static atoms and static versions of all the tree
* index matrices. The algorithm and algebra are documented in
* dynare++-ramsey.pdf. */
class PlannerBuilder {
friend class MultInitSS;
public:
/** Type for a set of variable names. */
typedef unordered_set<const char*> Tvarset;
/** Type for a set of equations. An equation is identified by
* an index to an equation in the equation vector given by
* DynareModel::eqs. The tree index of the i-th formula is
* retrieved as DynareModel::egs.formula(i). */
typedef vector<int> Teqset;
protected:
/** This is a set of variables wrt which the planner
* optimizes. These could be all endogenous variables, but it
* is beneficial to exclude all variables which are
* deterministic transformations of past exogenous variables,
* since the planner cannot influence them. This could save a
* few equations. This is not changed after it is constructed,
* but it is constructed manually, so it cannot be declared as
* const. */
Tvarset yset;
/** These are the equation indices constituing the constraints
* for the planner. Again, it is beneficial to exclude all
* equations defining exogenous variables excluded from
* yset. */
const Teqset fset;
/** Reference to the model. */
ogdyn::DynareModel& model;
/** Tree index of the planner objective. */
int tb;
/** Tree index of the planner discount parameter. */
int tbeta;
/** The maximum lead in the model including the planner's
* objective before building the planner's FOCs. */
const int maxlead;
/** The minimum lag in the model including the planner's objective
* before building the planner's FOCs. */
const int minlag;
/** Tree indices of formulas in the planner FOCs involving
* derivatives of the planner's objective. Rows correspond to the
* endogenous variables, columns correspond to lags in the
* objective function. The contents of the matrix will evolve as
* the algorithm proceeds. */
IntegerMatrix diff_b;
/** Tree indices of formulas in the planner FOCs involving
* derivatives of the model equations (constraints). The first
* dimension corresponds to endogenous variables, the second to
* the constraints, the third to lags or leads of endogenous
* variables in the constraints. The contents of the array will
* evolve as the algorithm proceeds.*/
IntegerArray3 diff_f;
/** Static version of the model atoms. It is needed to build
* static version of diff_b and diff_f. */
ogp::StaticFineAtoms static_atoms;
/** Static version of all the trees of diff_b and diff_f build
* over static_atoms. */
ogp::OperationTree static_tree;
/** Tree indices of static version of diff_b over static_atoms and static_tree. */
IntegerMatrix diff_b_static;
/** Tree indices of static version of diff_f over static_atoms
* and static_tree. This member is created before calling
* lagrange_mult_f(), so it does not contain the
* multiplication with the lagrange multipliers. */
IntegerArray3 diff_f_static;
/** Auxiliary variables mapping. During the algorithm, some
* auxiliary variables for the terms might be created, so we
* remember their names and tree indices of the terms. This
* maps a name to the tree index of an expression equal to the
* auxiliary variable at time zero. The auxiliary variables
* names point to the dynamic atoms storage, tree inidices to
* the dynamic model tree. */
Tsubstmap aux_map;
/** Static version of aux_map. The names point to static_atoms
* storage, the tree indices to the static_tree. */
Tsubstmap static_aux_map;
/** Information about the number of various things. */
PlannerInfo info;
public:
/** Build the planner problem for the given model optimizing
* through the given endogenous variables with the given
* constraints. We allow for a selection of a subset of
* equations and variables in order to eliminate exogenous
* predetermined process which cannot be influenced by the
* social planner. */
PlannerBuilder(ogdyn::DynareModel& m, const Tvarset& yyset,
const Teqset& ffset);
/** Construct a copy of the builder with provided model, which
* is supposed to be the copy of the model in the builder. */
PlannerBuilder(const PlannerBuilder& pb, ogdyn::DynareModel& m);
/** Return the information. */
const PlannerInfo& get_info() const
{return info;}
protected:
/** Differentiate the planner objective wrt endogenous
* variables with different lags. */
void add_derivatives_of_b();
/** Differentiate the constraints wrt endogenous variables
* with different lags and leads. */
void add_derivatives_of_f();
/** Shift derivatives of diff_b. */
void shift_derivatives_of_b();
/** Shift derivatives of diff_ff. */
void shift_derivatives_of_f();
/** Multiply with the discount factor terms in diff_b. */
void beta_multiply_b();
/** Multiply with the discount factor terms in diff_f. */
void beta_multiply_f();
/** Fill static_atoms and static_tree and build diff_b_static,
* diff_f_static and aux_map_static with static versions of diff_b,
* diff_f and aux_map. */
void make_static_version();
/** Multiply diff_f with Langrange multipliers. */
void lagrange_mult_f();
/** Add the equations to the mode, including equation for auxiliary variables. */
void form_equations();
private:
/** Fill yset for a given yyset and given name storage. */
void fill_yset(const ogp::NameStorage& ns, const Tvarset& yyset);
/** Fill aux_map and aux_map_static for a given aaux_map and
* aaux_map_static for a given storage of dynamic atoms (used
* for aux_map) and static atoms storage from this object for
* aux_map_static. */
void fill_aux_map(const ogp::NameStorage& ns, const Tsubstmap& aaux_map,
const Tsubstmap& astatic_aux_map);
/** Avoid copying from only PlannerBuilder. */
PlannerBuilder(const PlannerBuilder& pb);
};
/** This class only calculates for the given initial guess of
* endogenous variables, initial guess of the Langrange
* multipliers of the social planner problem yielding the least
* square error. It is used by just calling its constructor. The
* constructor takes non-const reference to the vector of
* endogenous variables, calculates lambdas and put the values of
* lambdas to the vector. The algbera is found in
* dynare++-ramsey.pdf.
*
* The code can be run only after the parsing has been finished in
* atoms. */
class MultInitSS : public ogp::FormulaEvalLoader {
protected:
/** The constant reference to the builder. */
const PlannerBuilder& builder;
/** The constant term of the problem. Its length is the number
* of endogenous variable wrt the planner optimizes. */
Vector b;
/** The matrix of the overdetermined problem. The number of
* rows is equal to the number of endogenous variables wrt
* which the planner optimizes, the number of columns is equal
* to the number of Langrange multipliers which is equal to
* the number of constraints which is smaller than the number
* of endogenous variables. Hence the system b+F*lambda=0 is
* overdetermined. */
GeneralMatrix F;
public:
/** The constructor of the object which does everything. Its
* main goal is to update yy. Note that if an item of yy
* corresponding to a lagrange multiplier is already set, it
* is not reset. */
MultInitSS(const PlannerBuilder& pb, const Vector& pvals, Vector& yy);
/** This loads evaluated parts of b or F and decodes i and
* advances b or F depending on the decoded i. The decoding is
* dependent on the way how the terms of builder.diff_b and
* builder.diff_f_save have been put the the
* ogp::FormulaCustomEvaluator. This is documented in the code
* of the constructor. */
void load(int i, double res);
};
};
#endif
// Local Variables:
// mode:C++
// End:
dynare++/sylv/Makefile.am deleted 100644 → 0
View file @ 0691b083
SUBDIRS = cc testing
EXTRA_DIST = sylvester.tex change_log.html matlab
if HAVE_PDFETEX
pdf-local: sylvester.pdf
sylvester.pdf: sylvester.tex
$(PDFETEX) sylvester
endif
CLEANFILES = *.pdf *.log
dynare++/sylv/cc/BlockDiagonal.cpp deleted 100644 → 0
View file @ 0691b083
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/BlockDiagonal.cpp,v 1.1.1.1 2004/06/04 13:00:20 kamenik Exp $ */
/* Tag $Name: $ */
#include "BlockDiagonal.h"
#include <cstdio>
#include <cstring>
BlockDiagonal::BlockDiagonal(const double* d, int d_size)
: QuasiTriangular(d, d_size),
row_len(new int[d_size]), col_len(new int[d_size])
{
for (int i = 0; i < d_size; i++) {
row_len[i] = d_size;
col_len[i] = 0;
}
}
BlockDiagonal::BlockDiagonal(const QuasiTriangular& t)
: QuasiTriangular(t),
row_len(new int[t.numRows()]), col_len(new int[t.numRows()])
{
for (int i = 0; i < t.numRows(); i++) {
row_len[i] = t.numRows();
col_len[i] = 0;
}
}
BlockDiagonal::BlockDiagonal(int p, const BlockDiagonal& b)
: QuasiTriangular(p, b),
row_len(new int[b.numRows()]), col_len(new int[b.numRows()])
{
memcpy(row_len, b.row_len, b.numRows()*sizeof(int));
memcpy(col_len, b.col_len, b.numRows()*sizeof(int));
}
BlockDiagonal::BlockDiagonal(const BlockDiagonal& b)
: QuasiTriangular(b),
row_len(new int[b.numRows()]), col_len(new int[b.numRows()])
{
memcpy(row_len, b.row_len, b.numRows()*sizeof(int));
memcpy(col_len, b.col_len, b.numRows()*sizeof(int));
}
/* put zeroes to right upper submatrix whose first column is defined
* by 'edge' */
void BlockDiagonal::setZerosToRU(diag_iter edge)
{
int iedge = (*edge).getIndex();
for (int i = 0; i < iedge; i++)
for (int j = iedge; j < numCols(); j++)
get(i,j) = 0.0;
}
/* Updates row_len and col_len so that there are zeroes in upper right part, this
* |T1 0 |
* |0 T2|. The first column of T2 is given by diagonal iterator 'edge'.
* Note the semantics of row_len and col_len. row_len[i] is distance
* of the right-most non-zero element of i-th row from the left, and
* col_len[j] is distance of top-most non-zero element of j-th column
* to the top. (First element has distance 1).
*/
void BlockDiagonal::setZeroBlockEdge(diag_iter edge)
{
setZerosToRU(edge);
int iedge = (*edge).getIndex();
for (diag_iter run = diag_begin(); run != edge; ++run) {
int ind = (*run).getIndex();
if (row_len[ind] > iedge) {
row_len[ind] = iedge;
if (!(*run).isReal())
row_len[ind+1] = iedge;
}
}
for (diag_iter run = edge; run != diag_end(); ++run) {
int ind = (*run).getIndex();
if (col_len[ind] < iedge) {
col_len[ind] = iedge;
if (!(*run).isReal())
col_len[ind+1] = iedge;
}
}
}
BlockDiagonal::const_col_iter
BlockDiagonal::col_begin(const DiagonalBlock& b) const
{
int jbar = b.getIndex();
int d_size = diagonal.getSize();
return const_col_iter(&getData()[jbar*d_size + col_len[jbar]], d_size,
b.isReal(), col_len[jbar]);
}
BlockDiagonal::col_iter
BlockDiagonal::col_begin(const DiagonalBlock& b)
{
int jbar = b.getIndex();
int d_size = diagonal.getSize();
return col_iter(&getData()[jbar*d_size + col_len[jbar]], d_size,
b.isReal(), col_len[jbar]);
}
BlockDiagonal::const_row_iter
BlockDiagonal::row_end(const DiagonalBlock& b) const
{
int jbar = b.getIndex();
int d_size = diagonal.getSize();
return const_row_iter(&getData()[d_size*row_len[jbar]+jbar], d_size,
b.isReal(), row_len[jbar]);
}
BlockDiagonal::row_iter
BlockDiagonal::row_end(const DiagonalBlock& b)
{
int jbar = b.getIndex();
int d_size = diagonal.getSize();
return row_iter(&getData()[d_size*row_len[jbar]+jbar], d_size,
b.isReal(), row_len[jbar]);
}
int BlockDiagonal::getNumZeros() const
{
int sum = 0;
for (int i = 0; i < diagonal.getSize(); i++) {
sum += diagonal.getSize() - row_len[i];
}
return sum;
}
QuasiTriangular::const_diag_iter
BlockDiagonal::findBlockStart(const_diag_iter from) const
{
if (from != diag_end()) {
++from;
while (from != diag_end() &&
col_len[(*from).getIndex()] != (*from).getIndex())
++from;
}
return from;
}
int BlockDiagonal::getLargestBlock() const
{
int largest = 0;
const_diag_iter start = diag_begin();
const_diag_iter end = findBlockStart(start);
while (start != diag_end()) {
int si = (*start).getIndex();
int ei = diagonal.getSize();
if (end != diag_end())
ei = (*end).getIndex();
if (largest < ei-si)
largest = ei-si;
start = end;
end = findBlockStart(start);
}
return largest;
}
void BlockDiagonal::savePartOfX(int si, int ei, const KronVector& x, Vector& work)
{
for (int i = si; i < ei; i++) {
ConstKronVector xi(x, i);
Vector target(work, (i-si)*xi.length(), xi.length());
target = xi;
}
}
void BlockDiagonal::multKronBlock(const_diag_iter start, const_diag_iter end,
KronVector& x, Vector& work) const
{
int si = (*start).getIndex();
int ei = diagonal.getSize();
if (end != diag_end())
ei = (*end).getIndex();
savePartOfX(si, ei, x, work);
for (const_diag_iter di = start; di != end; ++di) {
int jbar = (*di).getIndex();
if ((*di).isReal()) {
KronVector xi(x, jbar);
xi.zeros();
Vector wi(work, (jbar-si)*xi.length(), xi.length());
xi.add(*((*di).getAlpha()), wi);
for (const_row_iter ri = row_begin(*di); ri != row_end(*di); ++ri) {
int col = ri.getCol();
Vector wj(work, (col-si)*xi.length(), xi.length());
xi.add(*ri, wj);
}
} else {
KronVector xi(x, jbar);
KronVector xii(x, jbar+1);
xi.zeros();
xii.zeros();
Vector wi(work, (jbar-si)*xi.length(), xi.length());
Vector wii(work, (jbar+1-si)*xi.length(), xi.length());
xi.add(*((*di).getAlpha()), wi);
xi.add((*di).getBeta1(), wii);
xii.add((*di).getBeta2(), wi);
xii.add(*((*di).getAlpha()), wii);
for (const_row_iter ri = row_begin(*di); ri != row_end(*di); ++ri) {
int col = ri.getCol();
Vector wj(work, (col-si)*xi.length(), xi.length());
xi.add(ri.a(), wj);
xii.add(ri.b(), wj);
}
}
}
}
void BlockDiagonal::multKronBlockTrans(const_diag_iter start, const_diag_iter end,
KronVector& x, Vector& work) const
{
int si = (*start).getIndex();
int ei = diagonal.getSize();
if (end != diag_end())
ei = (*end).getIndex();
savePartOfX(si, ei, x, work);
for (const_diag_iter di = start; di != end; ++di) {
int jbar = (*di).getIndex();
if ((*di).isReal()) {
KronVector xi(x, jbar);
xi.zeros();
Vector wi(work, (jbar-si)*xi.length(), xi.length());
xi.add(*((*di).getAlpha()), wi);
for (const_col_iter ci = col_begin(*di); ci != col_end(*di); ++ci) {
int row = ci.getRow();
Vector wj(work, (row-si)*xi.length(), xi.length());
xi.add(*ci, wj);
}
} else {
KronVector xi(x, jbar);
KronVector xii(x, jbar+1);
xi.zeros();
xii.zeros();
Vector wi(work, (jbar-si)*xi.length(), xi.length());
Vector wii(work, (jbar+1-si)*xi.length(), xi.length());
xi.add(*((*di).getAlpha()), wi);
xi.add((*di).getBeta2(), wii);
xii.add((*di).getBeta1(), wi);
xii.add(*((*di).getAlpha()), wii);
for (const_col_iter ci = col_begin(*di); ci != col_end(*di); ++ci) {
int row = ci.getRow();
Vector wj(work, (row-si)*xi.length(), xi.length());
xi.add(ci.a(), wj);
xii.add(ci.b(), wj);
}
}
}
}
void BlockDiagonal::multKron(KronVector& x) const
{
int largest = getLargestBlock();
Vector work(largest*x.getN()*power(x.getM(),x.getDepth()-1));
const_diag_iter start = diag_begin();
const_diag_iter end = findBlockStart(start);
while (start != diag_end()) {
multKronBlock(start, end, x, work);
start = end;
end = findBlockStart(start);
}
}
void BlockDiagonal::multKronTrans(KronVector& x) const
{
int largest = getLargestBlock();
Vector work(largest*x.getN()*power(x.getM(),x.getDepth()-1));
const_diag_iter start = diag_begin();
const_diag_iter end = findBlockStart(start);
while (start != diag_end()) {
multKronBlockTrans(start, end, x, work);
start = end;
end = findBlockStart(start);
}
}
void BlockDiagonal::printInfo() const
{
printf("Block sizes:");
int num_blocks = 0;
const_diag_iter start = diag_begin();
const_diag_iter end = findBlockStart(start);
while (start != diag_end()) {
int si = (*start).getIndex();
int ei = diagonal.getSize();
if (end != diag_end())
ei = (*end).getIndex();
printf(" %d", ei-si);
num_blocks++;
start = end;
end = findBlockStart(start);
}
printf("\nNum blocks: %d\n", num_blocks);
printf("There are %d zeros out of %d\n",
getNumZeros(), getNumOffdiagonal());
}
int BlockDiagonal::getNumBlocks() const
{
int num_blocks = 0;
const_diag_iter start = diag_begin();
const_diag_iter end = findBlockStart(start);
while (start != diag_end()) {
num_blocks++;
start = end;
end = findBlockStart(start);
}
return num_blocks;
}
// Local Variables:
// mode:C++
// End:
dynare++/sylv/cc/BlockDiagonal.h deleted 100644 → 0
View file @ 0691b083
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/BlockDiagonal.h,v 1.1.1.1 2004/06/04 13:00:20 kamenik Exp $ */
/* Tag $Name: $ */
#ifndef BLOCK_DIAGONAL_H
#define BLOCK_DIAGONAL_H
#include "QuasiTriangular.h"
class BlockDiagonal : public QuasiTriangular {
int* const row_len;
int* const col_len;
public:
BlockDiagonal(const double* d, int d_size);
BlockDiagonal(int p, const BlockDiagonal& b);
BlockDiagonal(const BlockDiagonal& b);
BlockDiagonal(const QuasiTriangular& t);
const BlockDiagonal& operator=(const QuasiTriangular& t)
{GeneralMatrix::operator=(t); return *this;}
const BlockDiagonal& operator=(const BlockDiagonal& b);
~BlockDiagonal() {delete [] row_len; delete [] col_len;}
void setZeroBlockEdge(diag_iter edge);
int getNumZeros() const;
int getNumBlocks() const;
int getLargestBlock() const;
void printInfo() const;
void multKron(KronVector& x) const;
void multKronTrans(KronVector& x) const;
const_col_iter col_begin(const DiagonalBlock& b) const;
col_iter col_begin(const DiagonalBlock& b);
const_row_iter row_end(const DiagonalBlock& b) const;
row_iter row_end(const DiagonalBlock& b);
QuasiTriangular* clone() const
{return new BlockDiagonal(*this);}
private:
void setZerosToRU(diag_iter edge);
const_diag_iter findBlockStart(const_diag_iter from) const;
static void savePartOfX(int si, int ei, const KronVector& x, Vector& work);
void multKronBlock(const_diag_iter start, const_diag_iter end,
KronVector& x, Vector& work) const;
void multKronBlockTrans(const_diag_iter start, const_diag_iter end,
KronVector& x, Vector& work) const;
};
#endif /* BLOCK_DIAGONAL_H */
// Local Variables:
// mode:C++
// End:
dynare++/sylv/cc/GeneralMatrix.cpp deleted 100644 → 0
View file @ 0691b083
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/GeneralMatrix.cpp,v 1.4 2004/11/24 20:41:59 kamenik Exp $ */
/* Tag $Name: $ */
#include "SylvException.h"
#include "GeneralMatrix.h"
#include <dynblas.h>
#include <dynlapack.h>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cmath>
#include <limits>
int GeneralMatrix::md_length = 32;
GeneralMatrix::GeneralMatrix(const GeneralMatrix& m)
: data(m.rows*m.cols), rows(m.rows), cols(m.cols), ld(m.rows)
{
copy(m);
}
GeneralMatrix::GeneralMatrix(const ConstGeneralMatrix& m)
: data(m.rows*m.cols), rows(m.rows), cols(m.cols), ld(m.rows)
{
copy(m);
}
GeneralMatrix::GeneralMatrix(const GeneralMatrix& m, const char* dummy)
: data(m.rows*m.cols), rows(m.cols), cols(m.rows), ld(m.cols)
{
for (int i = 0; i < m.rows; i++)
for (int j = 0; j < m.cols; j++)
get(j,i) = m.get(i,j);
}
GeneralMatrix::GeneralMatrix(const ConstGeneralMatrix& m, const char* dummy)
: data(m.rows*m.cols), rows(m.cols), cols(m.rows), ld(m.cols)
{
for (int i = 0; i < m.rows; i++)
for (int j = 0; j < m.cols; j++)
get(j,i) = m.get(i,j);
}
GeneralMatrix::GeneralMatrix(const GeneralMatrix& m, int i, int j, int nrows, int ncols)
: data(nrows*ncols), rows(nrows), cols(ncols), ld(nrows)
{
copy(m, i, j);
}
GeneralMatrix::GeneralMatrix(GeneralMatrix& m, int i, int j, int nrows, int ncols)
: data(m.base()+m.ld*j+i, m.ld*(ncols-1)+nrows), rows(nrows), cols(ncols), ld(m.ld)
{}
GeneralMatrix::GeneralMatrix(const GeneralMatrix& a, const GeneralMatrix& b)
: data(a.rows*b.cols), rows(a.rows), cols(b.cols), ld(a.rows)
{
gemm("N", a, "N", b, 1.0, 0.0);
}
GeneralMatrix::GeneralMatrix(const GeneralMatrix& a, const GeneralMatrix& b, const char* dum)
: data(a.rows*b.rows), rows(a.rows), cols(b.rows), ld(a.rows)
{
gemm("N", a, "T", b, 1.0, 0.0);
}
GeneralMatrix::GeneralMatrix(const GeneralMatrix& a, const char* dum, const GeneralMatrix& b)
: data(a.cols*b.cols), rows(a.cols), cols(b.cols), ld(a.cols)
{
gemm("T", a, "N", b, 1.0, 0.0);
}
GeneralMatrix::GeneralMatrix(const GeneralMatrix& a, const char* dum1,
const GeneralMatrix& b, const char* dum2)
: data(a.cols*b.rows), rows(a.cols), cols(b.rows), ld(a.cols)
{
gemm("T", a, "T", b, 1.0, 0.0);
}
GeneralMatrix::~GeneralMatrix()
{
}
void GeneralMatrix::place(const ConstGeneralMatrix& m, int i, int j)
{
if (i + m.numRows() > numRows() ||
j + m.numCols() > numCols())
throw SYLV_MES_EXCEPTION("Bad submatrix placement, matrix dimensions exceeded.");
GeneralMatrix tmpsub(*this, i, j, m.numRows(), m.numCols());
tmpsub.copy(m);
}
void GeneralMatrix::mult(const ConstGeneralMatrix& a, const ConstGeneralMatrix& b)
{
gemm("N", a, "N", b, 1.0, 0.0);
}
void GeneralMatrix::multAndAdd(const ConstGeneralMatrix& a, const ConstGeneralMatrix& b,
double mult)
{
gemm("N", a, "N", b, mult, 1.0);
}
void GeneralMatrix::multAndAdd(const ConstGeneralMatrix& a, const ConstGeneralMatrix& b,
const char* dum, double mult)
{
gemm("N", a, "T", b, mult, 1.0);
}
void GeneralMatrix::multAndAdd(const ConstGeneralMatrix& a, const char* dum,
const ConstGeneralMatrix& b, double mult)
{
gemm("T", a, "N", b, mult, 1.0);
}
void GeneralMatrix::multAndAdd(const ConstGeneralMatrix& a, const char* dum1,
const ConstGeneralMatrix& b, const char* dum2, double mult)
{
gemm("T", a, "T", b, mult, 1.0);
}
void GeneralMatrix::addOuter(const ConstVector& a, double mult)
{
if (numRows() != numCols())
throw SYLV_MES_EXCEPTION("Matrix is not square in GeneralMatrix::addOuter.");
if (numRows() != a.length())
throw SYLV_MES_EXCEPTION("Wrong length of a vector in GeneralMatrix::addOuter.");
// since BLAS dsyr (symmetric rank 1 update) assumes symmetricity, we do this manually
for (int i = 0; i < numRows(); i++)
for (int j = i; j < numRows(); j++) {
double s = mult*a[i]*a[j];
get(i,j) = get(i,j) + s;
if (i != j)
get(j,i) = get(j,i) + s;
}
}
void GeneralMatrix::multRight(const ConstGeneralMatrix& m)
{
gemm_partial_right("N", m, 1.0, 0.0);
}
void GeneralMatrix::multLeft(const ConstGeneralMatrix& m)
{
gemm_partial_left("N", m, 1.0, 0.0);
}
void GeneralMatrix::multRightTrans(const ConstGeneralMatrix& m)
{
gemm_partial_right("T", m, 1.0, 0.0);
}
void GeneralMatrix::multLeftTrans(const ConstGeneralMatrix& m)
{
gemm_partial_left("T", m, 1.0, 0.0);
}
// here we must be careful for ld
void GeneralMatrix::zeros()
{
if (ld == rows)
data.zeros();
else {
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
get(i,j) = 0;
}
}
void GeneralMatrix::unit()
{
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
if (i == j)
get(i,j) = 1.0;
else
get(i,j) = 0.0;
}
void GeneralMatrix::nans()
{
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
get(i,j) = std::numeric_limits<double>::quiet_NaN();
}
void GeneralMatrix::infs()
{
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
get(i,j) = std::numeric_limits<double>::infinity();
}
// here we must be careful for ld
void GeneralMatrix::mult(double a)
{
if (ld == rows)
data.mult(a);
else {
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
get(i,j) *= a;
}
}
// here we must be careful for ld
void GeneralMatrix::add(double a, const ConstGeneralMatrix& m)
{
if (m.numRows() != rows || m.numCols() != cols)
throw SYLV_MES_EXCEPTION("Matrix has different size in GeneralMatrix::add.");
if (ld == rows && m.ld == m.rows)
data.add(a, m.data);
else {
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
get(i,j) += a*m.get(i,j);
}
}
void GeneralMatrix::add(double a, const ConstGeneralMatrix& m, const char* dum)
{
if (m.numRows() != cols || m.numCols() != rows)
throw SYLV_MES_EXCEPTION("Matrix has different size in GeneralMatrix::add.");
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
get(i,j) += a*m.get(j,i);
}
void GeneralMatrix::copy(const ConstGeneralMatrix& m, int ioff, int joff)
{
for (int i = 0; i < rows; i++)
for (int j = 0; j < cols; j++)
get(i,j) = m.get(i+ioff,j+joff);
}
void GeneralMatrix::gemm(const char* transa, const ConstGeneralMatrix& a,
const char* transb, const ConstGeneralMatrix& b,
double alpha, double beta)
{
int opa_rows = a.numRows();
int opa_cols = a.numCols();
if (!strcmp(transa, "T")) {
opa_rows = a.numCols();
opa_cols = a.numRows();
}
int opb_rows = b.numRows();
int opb_cols = b.numCols();
if (!strcmp(transb, "T")) {
opb_rows = b.numCols();
opb_cols = b.numRows();
}
if (opa_rows != numRows() ||
opb_cols != numCols() ||
opa_cols != opb_rows) {
throw SYLV_MES_EXCEPTION("Wrong dimensions for matrix multiplication.");
}
blas_int m = opa_rows;
blas_int n = opb_cols;
blas_int k = opa_cols;
blas_int lda = a.ld;
blas_int ldb = b.ld;
blas_int ldc = ld;
if (lda > 0 && ldb > 0 && ldc > 0) {
dgemm(transa, transb, &m, &n, &k, &alpha, a.data.base(), &lda,
b.data.base(), &ldb, &beta, data.base(), &ldc);
} else if (numRows()*numCols() > 0) {
if (beta == 0.0)
zeros();
else
mult(beta);
}
}
void GeneralMatrix::gemm_partial_left(const char* trans, const ConstGeneralMatrix& m,
double alpha, double beta)
{
int icol;
for (icol = 0; icol + md_length < cols; icol += md_length) {
GeneralMatrix incopy((const GeneralMatrix&)*this, 0, icol, rows, md_length);
GeneralMatrix inplace((GeneralMatrix&)*this, 0, icol, rows, md_length);
inplace.gemm(trans, m, "N", ConstGeneralMatrix(incopy), alpha, beta);
}
if (cols > icol) {
GeneralMatrix incopy((const GeneralMatrix&)*this, 0, icol, rows, cols - icol);
GeneralMatrix inplace((GeneralMatrix&)*this, 0, icol, rows, cols - icol);
inplace.gemm(trans, m, "N", ConstGeneralMatrix(incopy), alpha, beta);
}
}
void GeneralMatrix::gemm_partial_right(const char* trans, const ConstGeneralMatrix& m,
double alpha, double beta)
{
int irow;
for (irow = 0; irow + md_length < rows; irow += md_length) {
GeneralMatrix incopy((const GeneralMatrix&)*this, irow, 0, md_length, cols);
GeneralMatrix inplace((GeneralMatrix&)*this, irow, 0, md_length, cols);
inplace.gemm("N", ConstGeneralMatrix(incopy), trans, m, alpha, beta);
}
if (rows > irow) {
GeneralMatrix incopy((const GeneralMatrix&)*this, irow, 0, rows - irow, cols);
GeneralMatrix inplace((GeneralMatrix&)*this, irow, 0, rows - irow, cols);
inplace.gemm("N", ConstGeneralMatrix(incopy), trans, m, alpha, beta);
}
}
ConstGeneralMatrix::ConstGeneralMatrix(const GeneralMatrix& m, int i, int j, int nrows, int ncols)
: data(m.getData(), j*m.getLD()+i, (ncols-1)*m.getLD()+nrows), rows(nrows), cols(ncols), ld(m.getLD())
{
// can check that the submatirx is fully in the matrix
}
ConstGeneralMatrix::ConstGeneralMatrix(const ConstGeneralMatrix& m, int i, int j, int nrows, int ncols)
: data(m.getData(), j*m.getLD()+i, (ncols-1)*m.getLD()+nrows), rows(nrows), cols(ncols), ld(m.getLD())
{
// can check that the submatirx is fully in the matrix
}
ConstGeneralMatrix::ConstGeneralMatrix(const GeneralMatrix& m)
: data(m.data), rows(m.rows), cols(m.cols), ld(m.ld) {}
double ConstGeneralMatrix::getNormInf() const
{
double norm = 0.0;
for (int i = 0; i < numRows(); i++) {
ConstVector rowi(data.base()+i, ld, cols);
double normi = rowi.getNorm1();
if (norm < normi)
norm = normi;
}
return norm;
}
double ConstGeneralMatrix::getNorm1() const
{
double norm = 0.0;
for (int j = 0; j < numCols(); j++) {
ConstVector colj(data.base()+ld*j, 1, rows);
double normj = colj.getNorm1();
if (norm < normj)
norm = normj;
}
return norm;
}
void ConstGeneralMatrix::multVec(double a, Vector& x, double b, const ConstVector& d) const
{
if (x.length() != rows || cols != d.length()) {
throw SYLV_MES_EXCEPTION("Wrong dimensions for vector multiply.");
}
if (rows > 0) {
blas_int mm = rows;
blas_int nn = cols;
double alpha = b;
blas_int lda = ld;
blas_int incx = d.skip();
double beta = a;
blas_int incy = x.skip();
dgemv("N", &mm, &nn, &alpha, data.base(), &lda, d.base(), &incx,
&beta, x.base(), &incy);
}
}
void ConstGeneralMatrix::multVecTrans(double a, Vector& x, double b,
const ConstVector& d) const
{
if (x.length() != cols || rows != d.length()) {
throw SYLV_MES_EXCEPTION("Wrong dimensions for vector multiply.");
}
if (rows > 0) {
blas_int mm = rows;
blas_int nn = cols;
double alpha = b;
blas_int lda = rows;
blas_int incx = d.skip();
double beta = a;
blas_int incy = x.skip();
dgemv("T", &mm, &nn, &alpha, data.base(), &lda, d.base(), &incx,
&beta, x.base(), &incy);
}
}
/* m = inv(this)*m */
void ConstGeneralMatrix::multInvLeft(const char* trans, int mrows, int mcols, int mld, double* d) const
{
if (rows != cols) {
throw SYLV_MES_EXCEPTION("The matrix is not square for inversion.");
}
if (cols != mrows) {
throw SYLV_MES_EXCEPTION("Wrong dimensions for matrix inverse mutliply.");
}
if (rows > 0) {
GeneralMatrix inv(*this);
lapack_int* ipiv = new lapack_int[rows];
lapack_int info;
lapack_int rows2 = rows, mcols2 = mcols, mld2 = mld;
dgetrf(&rows2, &rows2, inv.getData().base(), &rows2, ipiv, &info);
dgetrs(trans, &rows2, &mcols2, inv.base(), &rows2, ipiv, d,
&mld2, &info);
delete [] ipiv;
}
}
/* m = inv(this)*m */
void ConstGeneralMatrix::multInvLeft(GeneralMatrix& m) const
{
multInvLeft("N", m.numRows(), m.numCols(), m.getLD(), m.getData().base());
}
/* m = inv(this')*m */
void ConstGeneralMatrix::multInvLeftTrans(GeneralMatrix& m) const
{
multInvLeft("T", m.numRows(), m.numCols(), m.getLD(), m.getData().base());
}
/* d = inv(this)*d */
void ConstGeneralMatrix::multInvLeft(Vector& d) const
{
if (d.skip() != 1) {
throw SYLV_MES_EXCEPTION("Skip!=1 not implemented in ConstGeneralMatrix::multInvLeft(Vector&)");
}
multInvLeft("N", d.length(), 1, d.length(), d.base());
}
/* d = inv(this')*d */
void ConstGeneralMatrix::multInvLeftTrans(Vector& d) const
{
if (d.skip() != 1) {
throw SYLV_MES_EXCEPTION("Skip!=1 not implemented in ConstGeneralMatrix::multInvLeft(Vector&)");
}
multInvLeft("T", d.length(), 1, d.length(), d.base());
}
bool ConstGeneralMatrix::isFinite() const
{
for (int i = 0; i < numRows(); i++)
for (int j = 0; j < numCols(); j++)
if (! std::isfinite(get(i,j)))
return false;
return true;
}
bool ConstGeneralMatrix::isZero() const
{
for (int i = 0; i < numRows(); i++)
for (int j = 0; j < numCols(); j++)
if (get(i,j) != 0.0)
return false;
return true;
}
void ConstGeneralMatrix::print() const
{
printf("rows=%d, cols=%d\n",rows, cols);
for (int i = 0; i < rows; i++) {
printf("row %d:\n",i);
for (int j = 0; j < cols; j++) {
printf("%6.3g ",get(i,j));
}
printf("\n");
}
}
dynare++/sylv/cc/GeneralMatrix.h deleted 100644 → 0
View file @ 0691b083
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/GeneralMatrix.h,v 1.3 2004/11/24 20:41:59 kamenik Exp $ */
/* Tag $Name: $ */
#ifndef GENERAL_MATRIX_H
#define GENERAL_MATRIX_H
#include "Vector.h"
class GeneralMatrix;
class ConstGeneralMatrix {
friend class GeneralMatrix;
protected:
ConstVector data;
int rows;
int cols;
int ld;
public:
ConstGeneralMatrix(const double* d, int m, int n)
: data(d, m*n), rows(m), cols(n), ld(m) {}
ConstGeneralMatrix(const GeneralMatrix& m);
ConstGeneralMatrix(const GeneralMatrix& m, int i, int j, int nrows, int ncols);
ConstGeneralMatrix(const ConstGeneralMatrix& m, int i, int j, int nrows, int ncols);
virtual ~ConstGeneralMatrix() {}
const double& get(int i, int j) const
{return data[j*ld+i];}
int numRows() const {return rows;}
int numCols() const {return cols;}
int getLD() const {return ld;}
const double* base() const {return data.base();}
const ConstVector& getData() const {return data;}
double getNormInf() const;
double getNorm1() const;
/* x = scalar(a)*x + scalar(b)*this*d */
void multVec(double a, Vector& x, double b, const ConstVector& d) const;
/* x = scalar(a)*x + scalar(b)*this'*d */
void multVecTrans(double a, Vector& x, double b, const ConstVector& d) const;
/* x = x + this*d */
void multaVec(Vector& x, const ConstVector& d) const
{multVec(1.0, x, 1.0, d);}
/* x = x + this'*d */
void multaVecTrans(Vector& x, const ConstVector& d) const
{multVecTrans(1.0, x, 1.0, d);}
/* x = x - this*d */
void multsVec(Vector& x, const ConstVector& d) const
{multVec(1.0, x, -1.0, d);}
/* x = x - this'*d */
void multsVecTrans(Vector& x, const ConstVector& d) const
{multVecTrans(1.0, x, -1.0, d);}
/* m = inv(this)*m */
void multInvLeft(GeneralMatrix& m) const;
/* m = inv(this')*m */
void multInvLeftTrans(GeneralMatrix& m) const;
/* d = inv(this)*d */
void multInvLeft(Vector& d) const;
/* d = inv(this')*d */
void multInvLeftTrans(Vector& d) const;
bool isFinite() const;
/** Returns true of the matrix is exactly zero. */
bool isZero() const;
virtual void print() const;
protected:
void multInvLeft(const char* trans, int mrows, int mcols, int mld, double* d) const;
};
class GeneralMatrix {
friend class ConstGeneralMatrix;
protected:
Vector data;
int rows;
int cols;
int ld;
public:
GeneralMatrix(int m, int n)
: data(m*n), rows(m), cols(n), ld(m) {}
GeneralMatrix(const double* d, int m, int n)
: data(d, m*n), rows(m), cols(n), ld(m) {}
GeneralMatrix(double* d, int m, int n)
: data(d, m*n), rows(m), cols(n), ld(m) {}
GeneralMatrix(const GeneralMatrix& m);
GeneralMatrix(const ConstGeneralMatrix& m);
GeneralMatrix(const GeneralMatrix&m, const char* dummy); // transpose
GeneralMatrix(const ConstGeneralMatrix&m, const char* dummy); // transpose
GeneralMatrix(const GeneralMatrix& m, int i, int j, int nrows, int ncols);
GeneralMatrix(GeneralMatrix& m, int i, int j, int nrows, int ncols);
/* this = a*b */
GeneralMatrix(const GeneralMatrix& a, const GeneralMatrix& b);
/* this = a*b' */
GeneralMatrix(const GeneralMatrix& a, const GeneralMatrix& b, const char* dum);
/* this = a'*b */
GeneralMatrix(const GeneralMatrix& a, const char* dum, const GeneralMatrix& b);
/* this = a'*b */
GeneralMatrix(const GeneralMatrix& a, const char* dum1,
const GeneralMatrix& b, const char* dum2);
virtual ~GeneralMatrix();
const GeneralMatrix& operator=(const GeneralMatrix& m)
{data=m.data; rows=m.rows; cols=m.cols; ld=m.ld; return *this;}
const double& get(int i, int j) const
{return data[j*ld+i];}
double& get(int i, int j)
{return data[j*ld+i];}
int numRows() const {return rows;}
int numCols() const {return cols;}
int getLD() const {return ld;}
double* base() {return data.base();}
const double* base() const {return data.base();}
Vector& getData() {return data;}
const Vector& getData() const {return data;}
double getNormInf() const
{return ConstGeneralMatrix(*this).getNormInf();}
double getNorm1() const
{return ConstGeneralMatrix(*this).getNorm1();}
/* place matrix m to the position (i,j) */
void place(const ConstGeneralMatrix& m, int i, int j);
void place(const GeneralMatrix& m, int i, int j)
{place(ConstGeneralMatrix(m), i, j);}
/* this = a*b */
void mult(const ConstGeneralMatrix& a, const ConstGeneralMatrix& b);
void mult(const GeneralMatrix& a, const GeneralMatrix& b)
{mult(ConstGeneralMatrix(a), ConstGeneralMatrix(b));}
/* this = this + scalar*a*b */
void multAndAdd(const ConstGeneralMatrix& a, const ConstGeneralMatrix& b,
double mult=1.0);
void multAndAdd(const GeneralMatrix& a, const GeneralMatrix& b,
double mult=1.0)
{multAndAdd(ConstGeneralMatrix(a), ConstGeneralMatrix(b), mult);}
/* this = this + scalar*a*b' */
void multAndAdd(const ConstGeneralMatrix& a, const ConstGeneralMatrix& b,
const char* dum, double mult=1.0);
void multAndAdd(const GeneralMatrix& a, const GeneralMatrix& b,
const char* dum, double mult=1.0)
{multAndAdd(ConstGeneralMatrix(a), ConstGeneralMatrix(b), dum, mult);}
/* this = this + scalar*a'*b */
void multAndAdd(const ConstGeneralMatrix& a, const char* dum, const ConstGeneralMatrix& b,
double mult=1.0);
void multAndAdd(const GeneralMatrix& a, const char* dum, const GeneralMatrix& b,
double mult=1.0)
{multAndAdd(ConstGeneralMatrix(a), dum, ConstGeneralMatrix(b), mult);}
/* this = this + scalar*a'*b' */
void multAndAdd(const ConstGeneralMatrix& a, const char* dum1,
const ConstGeneralMatrix& b, const char* dum2, double mult=1.0);
void multAndAdd(const GeneralMatrix& a, const char* dum1,
const GeneralMatrix& b, const char* dum2, double mult=1.0)
{multAndAdd(ConstGeneralMatrix(a), dum1, ConstGeneralMatrix(b),dum2, mult);}
/* this = this + scalar*a*a' */
void addOuter(const ConstVector& a, double mult=1.0);
void addOuter(const Vector& a, double mult=1.0)
{addOuter(ConstVector(a), mult);}
/* this = this * m */
void multRight(const ConstGeneralMatrix& m);
void multRight(const GeneralMatrix& m)
{multRight(ConstGeneralMatrix(m));}
/* this = m * this */
void multLeft(const ConstGeneralMatrix& m);
void multLeft(const GeneralMatrix& m)
{multLeft(ConstGeneralMatrix(m));}
/* this = this * m' */
void multRightTrans(const ConstGeneralMatrix& m);
void multRightTrans(const GeneralMatrix& m)
{multRightTrans(ConstGeneralMatrix(m));}
/* this = m' * this */
void multLeftTrans(const ConstGeneralMatrix& m);
void multLeftTrans(const GeneralMatrix& m)
{multLeftTrans(ConstGeneralMatrix(m));}
/* x = scalar(a)*x + scalar(b)*this*d */
void multVec(double a, Vector& x, double b, const ConstVector& d) const
{ConstGeneralMatrix(*this).multVec(a, x, b, d);}
/* x = scalar(a)*x + scalar(b)*this'*d */
void multVecTrans(double a, Vector& x, double b, const ConstVector& d) const
{ConstGeneralMatrix(*this).multVecTrans(a, x, b, d);}
/* x = x + this*d */
void multaVec(Vector& x, const ConstVector& d) const
{ConstGeneralMatrix(*this).multaVec(x, d);}
/* x = x + this'*d */
void multaVecTrans(Vector& x, const ConstVector& d) const
{ConstGeneralMatrix(*this).multaVecTrans(x, d);}
/* x = x - this*d */
void multsVec(Vector& x, const ConstVector& d) const
{ConstGeneralMatrix(*this).multsVec(x, d);}
/* x = x - this'*d */
void multsVecTrans(Vector& x, const ConstVector& d) const
{ConstGeneralMatrix(*this).multsVecTrans(x, d);}
/* this = zero */
void zeros();
/** this = unit (on main diagonal) */
void unit();
/* this = NaN */
void nans();
/* this = Inf */
void infs();
/* this = scalar*this */
void mult(double a);
/* this = this + scalar*m */
void add(double a, const ConstGeneralMatrix& m);
void add(double a, const GeneralMatrix& m)
{add(a, ConstGeneralMatrix(m));}
/* this = this + scalar*m' */
void add(double a, const ConstGeneralMatrix& m, const char* dum);
void add(double a, const GeneralMatrix& m, const char* dum)
{add(a, ConstGeneralMatrix(m), dum);}
bool isFinite() const
{return (ConstGeneralMatrix(*this)).isFinite();}
bool isZero() const
{return (ConstGeneralMatrix(*this)).isZero();}
virtual void print() const
{ConstGeneralMatrix(*this).print();}
private:
void copy(const ConstGeneralMatrix& m, int ioff = 0, int joff = 0);
void copy(const GeneralMatrix& m, int ioff = 0, int joff = 0)
{copy(ConstGeneralMatrix(m), ioff, joff);}
void gemm(const char* transa, const ConstGeneralMatrix& a,
const char* transb, const ConstGeneralMatrix& b,
double alpha, double beta);
void gemm(const char* transa, const GeneralMatrix& a,
const char* transb, const GeneralMatrix& b,
double alpha, double beta)
{gemm(transa, ConstGeneralMatrix(a), transb, ConstGeneralMatrix(b),
alpha, beta);}
/* this = this * op(m) (without whole copy of this) */
void gemm_partial_right(const char* trans, const ConstGeneralMatrix& m,
double alpha, double beta);
void gemm_partial_right(const char* trans, const GeneralMatrix& m,
double alpha, double beta)
{gemm_partial_right(trans, ConstGeneralMatrix(m), alpha, beta);}
/* this = op(m) *this (without whole copy of this) */
void gemm_partial_left(const char* trans, const ConstGeneralMatrix& m,
double alpha, double beta);
void gemm_partial_left(const char* trans, const GeneralMatrix& m,
double alpha, double beta)
{gemm_partial_left(trans, ConstGeneralMatrix(m), alpha, beta);}
/* number of rows/columns for copy used in gemm_partial_* */
static int md_length;
};
#endif /* GENERAL_MATRIX_H */
// Local Variables:
// mode:C++
// End:
dynare++/sylv/cc/GeneralSylvester.cpp deleted 100644 → 0
View file @ 0691b083
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/GeneralSylvester.cpp,v 1.1.1.1 2004/06/04 13:00:20 kamenik Exp $ */
/* Tag $Name: $ */
#include "GeneralSylvester.h"
#include "SchurDecomp.h"
#include "SylvException.h"
#include "TriangularSylvester.h"
#include "IterativeSylvester.h"
#include <ctime>
GeneralSylvester::GeneralSylvester(int ord, int n, int m, int zero_cols,
const double* da, const double* db,
const double* dc, const double* dd,
const SylvParams& ps)
: pars(ps),
mem_driver(pars, 1, m, n, ord), order(ord), a(da, n),
b(db, n, n-zero_cols), c(dc, m), d(dd, n, power(m, order)),
solved(false)
{
init();
}
GeneralSylvester::GeneralSylvester(int ord, int n, int m, int zero_cols,
const double* da, const double* db,
const double* dc, double* dd,
const SylvParams& ps)
: pars(ps),
mem_driver(pars, 0, m, n, ord), order(ord), a(da, n),
b(db, n, n-zero_cols), c(dc, m), d(dd, n, power(m, order)),
solved(false)
{
init();
}
GeneralSylvester::GeneralSylvester(int ord, int n, int m, int zero_cols,
const double* da, const double* db,
const double* dc, const double* dd,
bool alloc_for_check)
: pars(alloc_for_check),
mem_driver(pars, 1, m, n, ord), order(ord), a(da, n),
b(db, n, n-zero_cols), c(dc, m), d(dd, n, power(m, order)),
solved(false)
{
init();
}
GeneralSylvester::GeneralSylvester(int ord, int n, int m, int zero_cols,
const double* da, const double* db,
const double* dc, double* dd,
bool alloc_for_check)
: pars(alloc_for_check),
mem_driver(pars, 0, m, n, ord), order(ord), a(da, n),
b(db, n, n-zero_cols), c(dc, m), d(dd, n, power(m, order)),
solved(false)
{
init();
}
void GeneralSylvester::init()
{
GeneralMatrix ainvb(b);
double rcond1;
double rcondinf;
a.multInvLeft2(ainvb, d, rcond1, rcondinf);
pars.rcondA1 = rcond1;
pars.rcondAI = rcondinf;
bdecomp = new SchurDecompZero(ainvb);
cdecomp = new SimilarityDecomp(c.getData().base(), c.numRows(), *(pars.bs_norm));
cdecomp->check(pars, c);
cdecomp->infoToPars(pars);
if (*(pars.method) == SylvParams::recurse)
sylv = new TriangularSylvester(*bdecomp, *cdecomp);
else
sylv = new IterativeSylvester(*bdecomp, *cdecomp);
}
void GeneralSylvester::solve()
{
if (solved)
throw SYLV_MES_EXCEPTION("Attempt to run solve() more than once.");
mem_driver.setStackMode(true);
clock_t start = clock();
// multiply d
d.multLeftITrans(bdecomp->getQ());
d.multRightKron(cdecomp->getQ(), order);
// convert to KronVector
KronVector dkron(d.getData(), getM(), getN(), order);
// solve
sylv->solve(pars, dkron);
// multiply d back
d.multLeftI(bdecomp->getQ());
d.multRightKron(cdecomp->getInvQ(), order);
clock_t end = clock();
pars.cpu_time = ((double)(end-start))/CLOCKS_PER_SEC;
mem_driver.setStackMode(false);
solved = true;
}
void GeneralSylvester::check(const double* ds)
{
if (!solved)
throw SYLV_MES_EXCEPTION("Cannot run check on system, which is not solved yet.");
mem_driver.setStackMode(true);
// calculate xcheck = AX+BXC^i-D
SylvMatrix dcheck(d.numRows(), d.numCols());
dcheck.multLeft(b.numRows()-b.numCols(), b, d);
dcheck.multRightKron(c, order);
dcheck.multAndAdd(a,d);
ConstVector dv(ds, d.numRows()*d.numCols());
dcheck.getData().add(-1.0, dv);
// calculate relative norms
pars.mat_err1 = dcheck.getNorm1()/d.getNorm1();
pars.mat_errI = dcheck.getNormInf()/d.getNormInf();
pars.mat_errF = dcheck.getData().getNorm()/d.getData().getNorm();
pars.vec_err1 = dcheck.getData().getNorm1()/d.getData().getNorm1();
pars.vec_errI = dcheck.getData().getMax()/d.getData().getMax();
mem_driver.setStackMode(false);
}
GeneralSylvester::~GeneralSylvester()
{
delete bdecomp;
delete cdecomp;
delete sylv;
}
// Local Variables:
// mode:C++
// End:
dynare++/sylv/cc/GeneralSylvester.h deleted 100644 → 0
View file @ 0691b083
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/GeneralSylvester.h,v 1.1.1.1 2004/06/04 13:00:20 kamenik Exp $ */
/* Tag $Name: $ */
#ifndef GENERAL_SYLVESTER_H
#define GENERAL_SYLVESTER_H
#include "SylvMatrix.h"
#include "SylvMemory.h"
#include "SimilarityDecomp.h"
#include "SylvesterSolver.h"
class GeneralSylvester {
SylvParams pars;
SylvMemoryDriver mem_driver;
int order;
const SqSylvMatrix a;
const SylvMatrix b;
const SqSylvMatrix c;
SylvMatrix d;
bool solved;
SchurDecompZero* bdecomp;
SimilarityDecomp* cdecomp;
SylvesterSolver* sylv;
public:
/* construct with my copy of d*/
GeneralSylvester(int ord, int n, int m, int zero_cols,
const double* da, const double* db,
const double* dc, const double* dd,
const SylvParams& ps);
GeneralSylvester(int ord, int n, int m, int zero_cols,
const double* da, const double* db,
const double* dc, const double* dd,
bool alloc_for_check = false);
/* construct with provided storage for d */
GeneralSylvester(int ord, int n, int m, int zero_cols,
const double* da, const double* db,
const double* dc, double* dd,
bool alloc_for_check = false);
GeneralSylvester(int ord, int n, int m, int zero_cols,
const double* da, const double* db,
const double* dc, double* dd,
const SylvParams& ps);
virtual ~GeneralSylvester();
int getM() const {return c.numRows();}
int getN() const {return a.numRows();}
const double* getResult() const {return d.base();}
const SylvParams& getParams() const {return pars;}
SylvParams& getParams() {return pars;}
void solve();
void check(const double* ds);
private:
void init();
};
#endif /* GENERAL_SYLVESTER_H */
// Local Variables:
// mode:C++
// End:
dynare++/sylv/cc/IterativeSylvester.cpp deleted 100644 → 0
View file @ 0691b083
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/IterativeSylvester.cpp,v 1.1.1.1 2004/06/04 13:00:20 kamenik Exp $ */
/* Tag $Name: $ */
#include "IterativeSylvester.h"
#include "KronUtils.h"
void IterativeSylvester::solve(SylvParams& pars, KronVector& x) const
{
int max_steps = *(pars.max_num_iter);
int steps = 1;
double max_norm = *(pars.convergence_tol);
double norm = performFirstStep(x);
QuasiTriangular* kpow = matrixK->clone();
QuasiTriangular* fpow = matrixF->clone();
while (steps < max_steps && norm > max_norm) {
kpow->multRight(SqSylvMatrix(*kpow)); // be careful to make copy
fpow->multRight(SqSylvMatrix(*fpow)); // also here
norm = performStep(*kpow, *fpow, x);
steps++;
}
delete fpow;
delete kpow;
pars.converged = (norm <= max_norm);
pars.iter_last_norm = norm;
pars.num_iter = steps;
}
double IterativeSylvester::performFirstStep(KronVector& x) const
{
KronVector xtmp((const KronVector&)x);
KronUtils::multKron(*matrixF, *matrixK, xtmp);
x.add(-1., xtmp);
double norm = xtmp.getMax();
return norm;
}
double IterativeSylvester::performStep(const QuasiTriangular& k, const QuasiTriangular& f,
KronVector& x)
{
KronVector xtmp((const KronVector&)x);
KronUtils::multKron(f, k, xtmp);
x.add(1.0, xtmp);
double norm = xtmp.getMax();
return norm;
}
// Local Variables:
// mode:C++
// End:
dynare++/sylv/cc/IterativeSylvester.h deleted 100644 → 0
View file @ 0691b083
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/IterativeSylvester.h,v 1.1.1.1 2004/06/04 13:00:20 kamenik Exp $ */
/* Tag $Name: $ */
#ifndef ITERATIVE_SYLVESTER_H
#define ITERATIVE_SYLVESTER_H
#include "SylvesterSolver.h"
#include "KronVector.h"
#include "QuasiTriangular.h"
#include "SimilarityDecomp.h"
class IterativeSylvester : public SylvesterSolver {
public:
IterativeSylvester(const QuasiTriangular& k, const QuasiTriangular& f)
: SylvesterSolver(k, f) {}
IterativeSylvester(const SchurDecompZero& kdecomp, const SchurDecomp& fdecomp)
: SylvesterSolver(kdecomp, fdecomp) {}
IterativeSylvester(const SchurDecompZero& kdecomp, const SimilarityDecomp& fdecomp)
: SylvesterSolver(kdecomp, fdecomp) {}
void solve(SylvParams& pars, KronVector& x) const;
private:
double performFirstStep(KronVector& x) const;
static double performStep(const QuasiTriangular& k, const QuasiTriangular& f,
KronVector& x);
};
#endif /* ITERATIVE_SYLVESTER_H */
// Local Variables:
// mode:C++
// End:
dynare++/sylv/cc/KronUtils.cpp deleted 100644 → 0
View file @ 0691b083
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/KronUtils.cpp,v 1.1.1.1 2004/06/04 13:00:31 kamenik Exp $ */
/* Tag $Name: $ */
#include "KronUtils.h"
void KronUtils::multAtLevel(int level, const QuasiTriangular& t,
KronVector& x)
{
if (0 < level && level < x.getDepth()) {
for (int i = 0; i < x.getM(); i++) {
KronVector xi(x, i);
multAtLevel(level, t, xi);
}
} else if (0 == level && 0 < x.getDepth()) {
GeneralMatrix tmp(x.base(), x.getN(), power(x.getM(),x.getDepth()));
t.multLeftOther(tmp);
} else if (0 == level && 0 == x.getDepth()) {
Vector b((const Vector&)x);
t.multVec(x,b);
} else { // 0 < level == depth
t.multKron(x);
}
}
void KronUtils::multAtLevelTrans(int level, const QuasiTriangular& t,
KronVector& x)
{
if (0 < level && level < x.getDepth()) {
for (int i = 0; i < x.getM(); i++) {
KronVector xi(x, i);
multAtLevelTrans(level, t, xi);
}
} else if (0 == level && 0 < x.getDepth()) {
GeneralMatrix tmp(x.base(), x.getN(), power(x.getM(),x.getDepth()));
t.multLeftOtherTrans(tmp);
} else if (level == 0 && 0 == x.getDepth()) {
Vector b((const Vector&)x);
t.multVecTrans(x,b);
} else { // 0 < level == depth
t.multKronTrans(x);
}
}
void KronUtils::multKron(const QuasiTriangular& f, const QuasiTriangular& k,
KronVector& x)
{
multAtLevel(0, k, x);
if (x.getDepth() > 0) {
for (int level = 1; level <= x.getDepth(); level++)
multAtLevelTrans(level, f, x);
}
}
dynare++/sylv/cc/KronUtils.h deleted 100644 → 0
View file @ 0691b083
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/KronUtils.h,v 1.1.1.1 2004/06/04 13:00:31 kamenik Exp $ */
/* Tag $Name: $ */
#ifndef KRON_UTILS_H
#define KRON_UTILS_H
#include "KronVector.h"
#include "QuasiTriangular.h"
class KronUtils {
public:
/* multiplies I_m\otimes..\I_m\otimes T\otimes I_m...I_m\otimes I_n
with given b and returns x. T must be (m,m), number of
\otimes is b.getDepth(), level is a number of I_m's between T
and I_n plus 1. If level=0, then we multiply
\I_m\otimes ..\otimes I_m\otimes T, T is (n,n) */
static void multAtLevel(int level, const QuasiTriangular& t,
KronVector& x);
static void multAtLevelTrans(int level, const QuasiTriangular& t,
KronVector& x);
/* multiplies x=(F'\otimes F'\otimes..\otimes K)x */
static void multKron(const QuasiTriangular& f, const QuasiTriangular& k,
KronVector& x);
};
#endif /* KRON_UTILS_H */
// Local Variables:
// mode:C++
// End: