Estimation DLL, matrix library:
* fixed various bugs * added utility for computing the infinite norm * changed the direction of mat::transpose(A, B) for consistency with GSLtime-shift
parent
0ef2f8c4bd
commit
547a912d5e
|
@ -116,7 +116,8 @@ DecisionRules::compute(const Matrix &jacobian, Matrix &g_y, Matrix &g_u) throw (
|
|||
for (size_t j = 0; j < n_fwrd; j++)
|
||||
mat::col_copy(A, n_back_mixed + zeta_fwrd_mixed[pi_fwrd[j]], n_static, n - n_static,
|
||||
E, n_back_mixed + pi_fwrd[j], 0);
|
||||
mat::negate(MatrixView(E, 0, 0, n - n_static, n_fwrd + n_back + 2*n_mixed)); // Here we take the opposite of some of the zeros initialized in the constructor, but it is not a problem
|
||||
MatrixView E_tmp(E, 0, 0, n - n_static, n_fwrd + n_back + 2*n_mixed);
|
||||
mat::negate(E_tmp); // Here we take the opposite of some of the zeros initialized in the constructor, but it is not a problem
|
||||
|
||||
// Perform the generalized Schur
|
||||
size_t sdim;
|
||||
|
|
|
@ -27,6 +27,7 @@
|
|||
#include <cstdlib>
|
||||
#include <cassert>
|
||||
#include <cstring>
|
||||
#include <cmath>
|
||||
|
||||
#include "Vector.hh"
|
||||
|
||||
|
@ -285,28 +286,29 @@ namespace mat
|
|||
std::swap(M(i,j), M(j,i));
|
||||
}
|
||||
|
||||
//! Computes M2 = M1' (even for rectangular matrices)
|
||||
//! Computes M1 = M2' (even for rectangular matrices)
|
||||
template<class Mat1, class Mat2>
|
||||
inline void
|
||||
transpose(Mat1 &M1, Mat2 &M2)
|
||||
transpose(Mat1 &M1, const Mat2 &M2)
|
||||
{
|
||||
assert(M1.getRows() == M2.getCols() && M1.getCols() == M2.getRows());
|
||||
for (size_t i = 0; i < M1.getRows(); i++)
|
||||
for (size_t j = 0; j < M1.getCols(); j++)
|
||||
M2(j, i) = M1(i, j);
|
||||
M1(i, j) = M2(j, i);
|
||||
}
|
||||
|
||||
//! Computes m1 = m1 + m2
|
||||
template<class Mat1, class Mat2>
|
||||
void
|
||||
add(Mat1 m1, const Mat2 m2)
|
||||
add(Mat1 &m1, const Mat2 &m2)
|
||||
{
|
||||
assert(m1.getRows() == m2.getRows() && m1.getCols() == m2.getCols());
|
||||
double *p1 = m1.getData();
|
||||
const double *p2 = m2.getData();
|
||||
while (p1 < m1.getData() + m1.getCols() * m1.getLd())
|
||||
{
|
||||
double *pp1 = p1, *pp2 = p2;
|
||||
double *pp1 = p1;
|
||||
const double *pp2 = p2;
|
||||
while (pp1 < p1 + m1.getRows())
|
||||
*pp1++ += *pp2++;
|
||||
|
||||
|
@ -315,17 +317,18 @@ namespace mat
|
|||
}
|
||||
}
|
||||
|
||||
//! Computes m1 = m1 m2
|
||||
//! Computes m1 = m1 - m2
|
||||
template<class Mat1, class Mat2>
|
||||
void
|
||||
sub(Mat1 m1, const Mat2 m2)
|
||||
sub(Mat1 &m1, const Mat2 &m2)
|
||||
{
|
||||
assert(m1.getRows() == m2.getRows() && m1.getCols() == m2.getCols());
|
||||
double *p1 = m1.getData();
|
||||
const double *p2 = m2.getData();
|
||||
while (p1 < m1.getData() + m1.getCols() * m1.getLd())
|
||||
{
|
||||
double *pp1 = p1, *pp2 = p2;
|
||||
double *pp1 = p1;
|
||||
const double *pp2 = p2;
|
||||
while (pp1 < p1 + m1.getRows())
|
||||
*pp1++ -= *pp2++;
|
||||
|
||||
|
@ -337,7 +340,7 @@ namespace mat
|
|||
//! Does m = -m
|
||||
template<class Mat>
|
||||
void
|
||||
negate(Mat m)
|
||||
negate(Mat &m)
|
||||
{
|
||||
double *p = m.getData();
|
||||
while (p < m.getData() + m.getCols() * m.getLd())
|
||||
|
@ -353,6 +356,27 @@ namespace mat
|
|||
}
|
||||
}
|
||||
|
||||
// Computes the infinite norm of a matrix
|
||||
template<class Mat>
|
||||
double
|
||||
nrminf(const Mat &m)
|
||||
{
|
||||
double nrm = 0;
|
||||
const double *p = m.getData();
|
||||
while (p < m.getData() + m.getCols() * m.getLd())
|
||||
{
|
||||
const double *pp = p;
|
||||
while (pp < p + m.getRows())
|
||||
{
|
||||
if (fabs(*pp) > nrm)
|
||||
nrm = fabs(*pp);
|
||||
pp++;
|
||||
}
|
||||
|
||||
p += m.getLd();
|
||||
}
|
||||
return nrm;
|
||||
}
|
||||
} // End of namespace
|
||||
|
||||
#endif
|
||||
|
|
|
@ -25,6 +25,7 @@
|
|||
|
||||
#include <cstdlib>
|
||||
#include <cassert>
|
||||
#include <cmath>
|
||||
|
||||
/*
|
||||
This header defines three vector classes, which implement a "vector concept"
|
||||
|
@ -159,7 +160,7 @@ namespace vec
|
|||
//! Computes v1 = v1 + v2
|
||||
template<class Vec1, class Vec2>
|
||||
void
|
||||
add(Vec1 v1, const Vec2 v2)
|
||||
add(Vec1 &v1, const Vec2 &v2)
|
||||
{
|
||||
assert(v1.getSize() == v2.getSize());
|
||||
double *p1 = v1.getData();
|
||||
|
@ -175,7 +176,7 @@ namespace vec
|
|||
//! Computes v1 = v1 - v2
|
||||
template<class Vec1, class Vec2>
|
||||
void
|
||||
sub(Vec1 v1, const Vec2 v2)
|
||||
sub(Vec1 &v1, const Vec2 &v2)
|
||||
{
|
||||
assert(v1.getSize() == v2.getSize());
|
||||
double *p1 = v1.getData();
|
||||
|
@ -191,7 +192,7 @@ namespace vec
|
|||
//! Does v = -v
|
||||
template<class Vec>
|
||||
void
|
||||
negate(Vec v)
|
||||
negate(Vec &v)
|
||||
{
|
||||
double *p = v.getData();
|
||||
while (p < v.getData() + v.getSize() * v.getStride())
|
||||
|
@ -201,6 +202,21 @@ namespace vec
|
|||
}
|
||||
}
|
||||
|
||||
// Computes the infinite norm of a vector
|
||||
template<class Vec>
|
||||
double
|
||||
nrminf(const Vec &v)
|
||||
{
|
||||
double nrm = 0;
|
||||
const double *p = v.getData();
|
||||
while (p < v.getData() + v.getSize() * v.getStride())
|
||||
{
|
||||
if (fabs(*p) > nrm)
|
||||
nrm = fabs(*p);
|
||||
p += v.getStride();
|
||||
}
|
||||
return nrm;
|
||||
}
|
||||
} // End of namespace
|
||||
|
||||
#endif
|
||||
|
|
Loading…
Reference in New Issue