1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
|
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/SchurDecomp.cpp,v 1.1.1.1 2004/06/04 13:00:44 kamenik Exp $ */
/* Tag $Name: $ */
#include "SchurDecomp.h"
#include <dynlapack.h>
SchurDecomp::SchurDecomp(const SqSylvMatrix& m)
: q_destroy(true), t_destroy(true)
{
lapack_int rows = m.numRows();
q = new SqSylvMatrix(rows);
SqSylvMatrix auxt(m);
lapack_int sdim;
double* const wr = new double[rows];
double* const wi = new double[rows];
lapack_int lwork = 6*rows;
double* const work = new double[lwork];
lapack_int info;
dgees("V", "N", 0, &rows, auxt.base(), &rows, &sdim,
wr, wi, q->base(), &rows,
work, &lwork, 0, &info);
delete [] work;
delete [] wi;
delete [] wr;
t = new QuasiTriangular(auxt.base(), rows);
}
SchurDecomp::SchurDecomp(const QuasiTriangular& tr)
: q_destroy(true), t_destroy(true)
{
q = new SqSylvMatrix(tr.numRows());
q->setUnit();
t = new QuasiTriangular(tr);
}
SchurDecomp::SchurDecomp(QuasiTriangular& tr)
: q_destroy(true), t_destroy(false)
{
q = new SqSylvMatrix(tr.numRows());
q->setUnit();
t = &tr;
}
SchurDecomp::~SchurDecomp()
{
if (t_destroy)
delete t;
if (q_destroy)
delete q;
}
int SchurDecomp::getDim() const
{
return t->numRows();
}
SchurDecompZero::SchurDecompZero(const GeneralMatrix& m)
: SchurDecomp(SqSylvMatrix(m, m.numRows()-m.numCols(), 0, m.numCols())),
ru(m, 0, 0, m.numRows()-m.numCols(), m.numCols())
{
ru.multRight(getQ());
}
int SchurDecompZero::getDim() const
{
return getT().numRows()+ru.numRows();
}
|
