1 /*
2 * Copyright © 2004-2011 Ondra Kamenik
3 * Copyright © 2019 Dynare Team
4 *
5 * This file is part of Dynare.
6 *
7 * Dynare is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * Dynare is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with Dynare. If not, see <http://www.gnu.org/licenses/>.
19 */
20
21 #include "GeneralSylvester.hh"
22 #include "SchurDecomp.hh"
23 #include "SylvException.hh"
24 #include "TriangularSylvester.hh"
25 #include "IterativeSylvester.hh"
26 #include "int_power.hh"
27
28 #include <ctime>
29
GeneralSylvester(int ord,int n,int m,int zero_cols,const ConstVector & da,const ConstVector & db,const ConstVector & dc,const ConstVector & dd,const SylvParams & ps)30 GeneralSylvester::GeneralSylvester(int ord, int n, int m, int zero_cols,
31 const ConstVector &da, const ConstVector &db,
32 const ConstVector &dc, const ConstVector &dd,
33 const SylvParams &ps)
34 : pars(ps),
35 order(ord), a(Vector{da}, n),
36 b(Vector{db}, n, n-zero_cols), c(Vector{dc}, m), d(Vector{dd}, n, power(m, order)),
37 solved(false)
38 {
39 init();
40 }
41
GeneralSylvester(int ord,int n,int m,int zero_cols,const ConstVector & da,const ConstVector & db,const ConstVector & dc,Vector & dd,const SylvParams & ps)42 GeneralSylvester::GeneralSylvester(int ord, int n, int m, int zero_cols,
43 const ConstVector &da, const ConstVector &db,
44 const ConstVector &dc, Vector &dd,
45 const SylvParams &ps)
46 : pars(ps),
47 order(ord), a(Vector{da}, n),
48 b(Vector{db}, n, n-zero_cols), c(Vector{dc}, m), d(dd, n, power(m, order)),
49 solved(false)
50 {
51 init();
52 }
53
GeneralSylvester(int ord,int n,int m,int zero_cols,const ConstVector & da,const ConstVector & db,const ConstVector & dc,const ConstVector & dd,bool alloc_for_check)54 GeneralSylvester::GeneralSylvester(int ord, int n, int m, int zero_cols,
55 const ConstVector &da, const ConstVector &db,
56 const ConstVector &dc, const ConstVector &dd,
57 bool alloc_for_check)
58 : pars(alloc_for_check),
59 order(ord), a(Vector{da}, n),
60 b(Vector{db}, n, n-zero_cols), c(Vector{dc}, m), d(Vector{dd}, n, power(m, order)),
61 solved(false)
62 {
63 init();
64 }
65
GeneralSylvester(int ord,int n,int m,int zero_cols,const ConstVector & da,const ConstVector & db,const ConstVector & dc,Vector & dd,bool alloc_for_check)66 GeneralSylvester::GeneralSylvester(int ord, int n, int m, int zero_cols,
67 const ConstVector &da, const ConstVector &db,
68 const ConstVector &dc, Vector &dd,
69 bool alloc_for_check)
70 : pars(alloc_for_check),
71 order(ord), a(Vector{da}, n),
72 b(Vector{db}, n, n-zero_cols), c(Vector{dc}, m), d(dd, n, power(m, order)),
73 solved(false)
74 {
75 init();
76 }
77
78 void
init()79 GeneralSylvester::init()
80 {
81 GeneralMatrix ainvb(b);
82 double rcond1;
83 double rcondinf;
84 a.multInvLeft2(ainvb, d, rcond1, rcondinf);
85 pars.rcondA1 = rcond1;
86 pars.rcondAI = rcondinf;
87 bdecomp = std::make_unique<SchurDecompZero>(ainvb);
88 cdecomp = std::make_unique<SimilarityDecomp>(c.getData(), c.nrows(), *(pars.bs_norm));
89 cdecomp->check(pars, c);
90 cdecomp->infoToPars(pars);
91 if (*(pars.method) == SylvParams::solve_method::recurse)
92 sylv = std::make_unique<TriangularSylvester>(*bdecomp, *cdecomp);
93 else
94 sylv = std::make_unique<IterativeSylvester>(*bdecomp, *cdecomp);
95 }
96
97 void
solve()98 GeneralSylvester::solve()
99 {
100 if (solved)
101 throw SYLV_MES_EXCEPTION("Attempt to run solve() more than once.");
102
103 clock_t start = clock();
104 // multiply d
105 d.multLeftITrans(bdecomp->getQ());
106 d.multRightKron(cdecomp->getQ(), order);
107 // convert to KronVector
108 KronVector dkron(d.getData(), getM(), getN(), order);
109 // solve
110 sylv->solve(pars, dkron);
111 // multiply d back
112 d.multLeftI(bdecomp->getQ());
113 d.multRightKron(cdecomp->getInvQ(), order);
114 clock_t end = clock();
115 pars.cpu_time = static_cast<double>(end-start)/CLOCKS_PER_SEC;
116
117 solved = true;
118 }
119
120 void
check(const ConstVector & ds)121 GeneralSylvester::check(const ConstVector &ds)
122 {
123 if (!solved)
124 throw SYLV_MES_EXCEPTION("Cannot run check on system, which is not solved yet.");
125
126 // calculate xcheck = A·X+B·X·⊗ⁱC−D
127 SylvMatrix dcheck(d.nrows(), d.ncols());
128 dcheck.multLeft(b.nrows()-b.ncols(), b, d);
129 dcheck.multRightKron(c, order);
130 dcheck.multAndAdd(a, d);
131 dcheck.getData().add(-1.0, ds);
132 // calculate relative norms
133 pars.mat_err1 = dcheck.getNorm1()/d.getNorm1();
134 pars.mat_errI = dcheck.getNormInf()/d.getNormInf();
135 pars.mat_errF = dcheck.getData().getNorm()/d.getData().getNorm();
136 pars.vec_err1 = dcheck.getData().getNorm1()/d.getData().getNorm1();
137 pars.vec_errI = dcheck.getData().getMax()/d.getData().getMax();
138 }
139