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
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
|
/* $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:
|
