2019-01-08 17:12:05 +01:00
|
|
|
/* $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.hh"
|
|
|
|
#include "GeneralMatrix.hh"
|
|
|
|
|
|
|
|
#include <dynblas.h>
|
|
|
|
#include <dynlapack.h>
|
|
|
|
|
2019-01-16 17:52:16 +01:00
|
|
|
#include <iostream>
|
|
|
|
#include <iomanip>
|
2019-01-08 17:12:05 +01:00
|
|
|
#include <cstdlib>
|
|
|
|
#include <cmath>
|
|
|
|
#include <limits>
|
2019-01-16 17:52:16 +01:00
|
|
|
#include <vector>
|
2019-01-08 17:12:05 +01:00
|
|
|
|
|
|
|
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, 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)
|
Dynare++ / sylvester equation solver: refactor Vector and ConstVector classes
- these classes now encapsulate a std::shared_ptr<{const, }double>, so that
they do not perform memory management, and several {Const,}Vector instances
can transparently share the same underlying data
- make converting constructor from ConstVector to Vector explicit, since that
entails memory allocation (but the reverse conversion is almost costless, so
keep it implicit); do the same for GeneralMatrix/ConstGeneralMatrix,
TwoDMatrix/ConstTwoDMatrix
- remove the constructors that were extracting a row/column from a matrix, and
replace them by getRow() and getCol() methods on {Const,}GeneralMatrix
- rename and change the API of the complex version Vector::add(), so that it is
explicit that it deals with complex numbers
- add constructors that take a MATLAB mxArray
2019-01-22 16:07:44 +01:00
|
|
|
: data(m.data, m.ld*j+i, m.ld*(ncols-1)+nrows), rows(nrows), cols(ncols), ld(m.ld)
|
2019-01-08 17:12:05 +01:00
|
|
|
{
|
|
|
|
}
|
|
|
|
|
2019-01-24 15:22:36 +01:00
|
|
|
GeneralMatrix &
|
|
|
|
GeneralMatrix::operator=(const ConstGeneralMatrix &m)
|
|
|
|
{
|
|
|
|
data = m.data;
|
|
|
|
rows = m.rows;
|
|
|
|
cols = m.cols;
|
|
|
|
ld = m.ld;
|
|
|
|
return *this;
|
|
|
|
}
|
|
|
|
|
Dynare++ / sylvester equation solver: refactor Vector and ConstVector classes
- these classes now encapsulate a std::shared_ptr<{const, }double>, so that
they do not perform memory management, and several {Const,}Vector instances
can transparently share the same underlying data
- make converting constructor from ConstVector to Vector explicit, since that
entails memory allocation (but the reverse conversion is almost costless, so
keep it implicit); do the same for GeneralMatrix/ConstGeneralMatrix,
TwoDMatrix/ConstTwoDMatrix
- remove the constructors that were extracting a row/column from a matrix, and
replace them by getRow() and getCol() methods on {Const,}GeneralMatrix
- rename and change the API of the complex version Vector::add(), so that it is
explicit that it deals with complex numbers
- add constructors that take a MATLAB mxArray
2019-01-22 16:07:44 +01:00
|
|
|
Vector
|
|
|
|
GeneralMatrix::getRow(int row)
|
|
|
|
{
|
|
|
|
return Vector{data, row, ld, cols};
|
|
|
|
}
|
|
|
|
|
|
|
|
Vector
|
|
|
|
GeneralMatrix::getCol(int col)
|
|
|
|
{
|
|
|
|
return Vector{data, ld*col, rows};
|
|
|
|
}
|
|
|
|
|
|
|
|
ConstVector
|
|
|
|
GeneralMatrix::getRow(int row) const
|
|
|
|
{
|
|
|
|
return ConstVector{data, row, ld, cols};
|
|
|
|
}
|
|
|
|
|
|
|
|
ConstVector
|
|
|
|
GeneralMatrix::getCol(int col) const
|
|
|
|
{
|
|
|
|
return ConstVector{data, ld*col, rows};
|
|
|
|
}
|
|
|
|
|
2019-01-08 17:12:05 +01:00
|
|
|
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,
|
2019-01-25 15:27:20 +01:00
|
|
|
const std::string &dum, double mult)
|
2019-01-08 17:12:05 +01:00
|
|
|
{
|
|
|
|
gemm("N", a, "T", b, mult, 1.0);
|
|
|
|
}
|
|
|
|
|
|
|
|
void
|
2019-01-25 15:27:20 +01:00
|
|
|
GeneralMatrix::multAndAdd(const ConstGeneralMatrix &a, const std::string &dum,
|
2019-01-08 17:12:05 +01:00
|
|
|
const ConstGeneralMatrix &b, double mult)
|
|
|
|
{
|
|
|
|
gemm("T", a, "N", b, mult, 1.0);
|
|
|
|
}
|
|
|
|
|
|
|
|
void
|
2019-01-25 15:27:20 +01:00
|
|
|
GeneralMatrix::multAndAdd(const ConstGeneralMatrix &a, const std::string &dum1,
|
|
|
|
const ConstGeneralMatrix &b, const std::string &dum2, double mult)
|
2019-01-08 17:12:05 +01:00
|
|
|
{
|
|
|
|
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
|
2019-01-16 17:52:16 +01:00
|
|
|
for (int i = 0; i < rows; i++)
|
|
|
|
for (int j = 0; j < cols; j++)
|
|
|
|
get(i, j) = 0;
|
2019-01-08 17:12:05 +01:00
|
|
|
}
|
|
|
|
|
|
|
|
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
|
2019-01-16 17:52:16 +01:00
|
|
|
for (int i = 0; i < rows; i++)
|
|
|
|
for (int j = 0; j < cols; j++)
|
|
|
|
get(i, j) *= a;
|
2019-01-08 17:12:05 +01:00
|
|
|
}
|
|
|
|
|
|
|
|
// 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
|
2019-01-16 17:52:16 +01:00
|
|
|
for (int i = 0; i < rows; i++)
|
|
|
|
for (int j = 0; j < cols; j++)
|
|
|
|
get(i, j) += a*m.get(i, j);
|
2019-01-08 17:12:05 +01:00
|
|
|
}
|
|
|
|
|
|
|
|
void
|
2019-01-25 15:27:20 +01:00
|
|
|
GeneralMatrix::add(double a, const ConstGeneralMatrix &m, const std::string &dum)
|
2019-01-08 17:12:05 +01:00
|
|
|
{
|
|
|
|
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
|
2019-01-25 15:27:20 +01:00
|
|
|
GeneralMatrix::gemm(const std::string &transa, const ConstGeneralMatrix &a,
|
|
|
|
const std::string &transb, const ConstGeneralMatrix &b,
|
2019-01-08 17:12:05 +01:00
|
|
|
double alpha, double beta)
|
|
|
|
{
|
|
|
|
int opa_rows = a.numRows();
|
|
|
|
int opa_cols = a.numCols();
|
2019-01-25 15:27:20 +01:00
|
|
|
if (transa == "T")
|
2019-01-08 17:12:05 +01:00
|
|
|
{
|
|
|
|
opa_rows = a.numCols();
|
|
|
|
opa_cols = a.numRows();
|
|
|
|
}
|
|
|
|
int opb_rows = b.numRows();
|
|
|
|
int opb_cols = b.numCols();
|
2019-01-25 15:27:20 +01:00
|
|
|
if (transb == "T")
|
2019-01-08 17:12:05 +01:00
|
|
|
{
|
|
|
|
opb_rows = b.numCols();
|
|
|
|
opb_cols = b.numRows();
|
|
|
|
}
|
|
|
|
|
|
|
|
if (opa_rows != numRows()
|
|
|
|
|| opb_cols != numCols()
|
|
|
|
|| opa_cols != opb_rows)
|
2019-03-08 15:32:13 +01:00
|
|
|
throw SYLV_MES_EXCEPTION("Wrong dimensions for matrix multiplication.");
|
2019-01-08 17:12:05 +01:00
|
|
|
|
|
|
|
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)
|
2019-03-08 15:32:13 +01:00
|
|
|
dgemm(transa.c_str(), transb.c_str(), &m, &n, &k, &alpha, a.data.base(), &lda,
|
|
|
|
b.data.base(), &ldb, &beta, data.base(), &ldc);
|
2019-01-08 17:12:05 +01:00
|
|
|
else if (numRows()*numCols() > 0)
|
|
|
|
{
|
|
|
|
if (beta == 0.0)
|
|
|
|
zeros();
|
|
|
|
else
|
|
|
|
mult(beta);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
void
|
2019-01-25 15:27:20 +01:00
|
|
|
GeneralMatrix::gemm_partial_left(const std::string &trans, const ConstGeneralMatrix &m,
|
2019-01-08 17:12:05 +01:00
|
|
|
double alpha, double beta)
|
|
|
|
{
|
|
|
|
int icol;
|
|
|
|
for (icol = 0; icol + md_length < cols; icol += md_length)
|
|
|
|
{
|
2019-03-08 15:32:13 +01:00
|
|
|
GeneralMatrix incopy(const_cast<const GeneralMatrix &>(*this), 0, icol, rows, md_length);
|
|
|
|
GeneralMatrix inplace(*this, 0, icol, rows, md_length);
|
2019-01-08 17:12:05 +01:00
|
|
|
inplace.gemm(trans, m, "N", ConstGeneralMatrix(incopy), alpha, beta);
|
|
|
|
}
|
|
|
|
if (cols > icol)
|
|
|
|
{
|
2019-03-08 15:32:13 +01:00
|
|
|
GeneralMatrix incopy(const_cast<const GeneralMatrix &>(*this), 0, icol, rows, cols - icol);
|
|
|
|
GeneralMatrix inplace(*this, 0, icol, rows, cols - icol);
|
2019-01-08 17:12:05 +01:00
|
|
|
inplace.gemm(trans, m, "N", ConstGeneralMatrix(incopy), alpha, beta);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
void
|
2019-01-25 15:27:20 +01:00
|
|
|
GeneralMatrix::gemm_partial_right(const std::string &trans, const ConstGeneralMatrix &m,
|
2019-01-08 17:12:05 +01:00
|
|
|
double alpha, double beta)
|
|
|
|
{
|
|
|
|
int irow;
|
|
|
|
for (irow = 0; irow + md_length < rows; irow += md_length)
|
|
|
|
{
|
2019-03-08 15:32:13 +01:00
|
|
|
GeneralMatrix incopy(const_cast<const GeneralMatrix &>(*this), irow, 0, md_length, cols);
|
|
|
|
GeneralMatrix inplace(*this, irow, 0, md_length, cols);
|
2019-01-08 17:12:05 +01:00
|
|
|
inplace.gemm("N", ConstGeneralMatrix(incopy), trans, m, alpha, beta);
|
|
|
|
}
|
|
|
|
if (rows > irow)
|
|
|
|
{
|
2019-03-08 15:32:13 +01:00
|
|
|
GeneralMatrix incopy(const_cast<const GeneralMatrix &>(*this), irow, 0, rows - irow, cols);
|
|
|
|
GeneralMatrix inplace(*this, irow, 0, rows - irow, cols);
|
2019-01-08 17:12:05 +01:00
|
|
|
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())
|
|
|
|
{
|
2019-03-27 19:22:35 +01:00
|
|
|
// FIXME: check that the submatrix is fully in the matrix
|
2019-01-08 17:12:05 +01:00
|
|
|
}
|
|
|
|
|
|
|
|
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())
|
|
|
|
{
|
2019-03-27 19:22:35 +01:00
|
|
|
// FIXME: check that the submatrix is fully in the matrix
|
2019-01-08 17:12:05 +01:00
|
|
|
}
|
|
|
|
|
|
|
|
ConstGeneralMatrix::ConstGeneralMatrix(const GeneralMatrix &m)
|
|
|
|
: data(m.data), rows(m.rows), cols(m.cols), ld(m.ld)
|
|
|
|
{
|
|
|
|
}
|
|
|
|
|
Dynare++ / sylvester equation solver: refactor Vector and ConstVector classes
- these classes now encapsulate a std::shared_ptr<{const, }double>, so that
they do not perform memory management, and several {Const,}Vector instances
can transparently share the same underlying data
- make converting constructor from ConstVector to Vector explicit, since that
entails memory allocation (but the reverse conversion is almost costless, so
keep it implicit); do the same for GeneralMatrix/ConstGeneralMatrix,
TwoDMatrix/ConstTwoDMatrix
- remove the constructors that were extracting a row/column from a matrix, and
replace them by getRow() and getCol() methods on {Const,}GeneralMatrix
- rename and change the API of the complex version Vector::add(), so that it is
explicit that it deals with complex numbers
- add constructors that take a MATLAB mxArray
2019-01-22 16:07:44 +01:00
|
|
|
ConstVector
|
|
|
|
ConstGeneralMatrix::getRow(int row) const
|
|
|
|
{
|
|
|
|
return ConstVector{data, row, ld, cols};
|
|
|
|
}
|
|
|
|
|
|
|
|
ConstVector
|
|
|
|
ConstGeneralMatrix::getCol(int col) const
|
|
|
|
{
|
|
|
|
return ConstVector{data, ld*col, rows};
|
|
|
|
}
|
|
|
|
|
2019-01-08 17:12:05 +01:00
|
|
|
double
|
|
|
|
ConstGeneralMatrix::getNormInf() const
|
|
|
|
{
|
|
|
|
double norm = 0.0;
|
|
|
|
for (int i = 0; i < numRows(); i++)
|
|
|
|
{
|
Dynare++ / sylvester equation solver: refactor Vector and ConstVector classes
- these classes now encapsulate a std::shared_ptr<{const, }double>, so that
they do not perform memory management, and several {Const,}Vector instances
can transparently share the same underlying data
- make converting constructor from ConstVector to Vector explicit, since that
entails memory allocation (but the reverse conversion is almost costless, so
keep it implicit); do the same for GeneralMatrix/ConstGeneralMatrix,
TwoDMatrix/ConstTwoDMatrix
- remove the constructors that were extracting a row/column from a matrix, and
replace them by getRow() and getCol() methods on {Const,}GeneralMatrix
- rename and change the API of the complex version Vector::add(), so that it is
explicit that it deals with complex numbers
- add constructors that take a MATLAB mxArray
2019-01-22 16:07:44 +01:00
|
|
|
double normi = getRow(i).getNorm1();
|
2019-03-08 15:32:13 +01:00
|
|
|
norm = std::max(normi, norm);
|
2019-01-08 17:12:05 +01:00
|
|
|
}
|
|
|
|
return norm;
|
|
|
|
}
|
|
|
|
|
|
|
|
double
|
|
|
|
ConstGeneralMatrix::getNorm1() const
|
|
|
|
{
|
|
|
|
double norm = 0.0;
|
|
|
|
for (int j = 0; j < numCols(); j++)
|
|
|
|
{
|
Dynare++ / sylvester equation solver: refactor Vector and ConstVector classes
- these classes now encapsulate a std::shared_ptr<{const, }double>, so that
they do not perform memory management, and several {Const,}Vector instances
can transparently share the same underlying data
- make converting constructor from ConstVector to Vector explicit, since that
entails memory allocation (but the reverse conversion is almost costless, so
keep it implicit); do the same for GeneralMatrix/ConstGeneralMatrix,
TwoDMatrix/ConstTwoDMatrix
- remove the constructors that were extracting a row/column from a matrix, and
replace them by getRow() and getCol() methods on {Const,}GeneralMatrix
- rename and change the API of the complex version Vector::add(), so that it is
explicit that it deals with complex numbers
- add constructors that take a MATLAB mxArray
2019-01-22 16:07:44 +01:00
|
|
|
double normj = getCol(j).getNorm1();
|
2019-03-08 15:32:13 +01:00
|
|
|
norm = std::max(normj, norm);
|
2019-01-08 17:12:05 +01:00
|
|
|
}
|
|
|
|
return norm;
|
|
|
|
}
|
|
|
|
|
|
|
|
void
|
|
|
|
ConstGeneralMatrix::multVec(double a, Vector &x, double b, const ConstVector &d) const
|
|
|
|
{
|
|
|
|
if (x.length() != rows || cols != d.length())
|
2019-01-16 17:52:16 +01:00
|
|
|
throw SYLV_MES_EXCEPTION("Wrong dimensions for vector multiply.");
|
2019-01-08 17:12:05 +01:00
|
|
|
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())
|
2019-01-16 17:52:16 +01:00
|
|
|
throw SYLV_MES_EXCEPTION("Wrong dimensions for vector multiply.");
|
2019-01-08 17:12:05 +01:00
|
|
|
if (rows > 0)
|
|
|
|
{
|
|
|
|
blas_int mm = rows;
|
|
|
|
blas_int nn = cols;
|
|
|
|
double alpha = b;
|
2019-01-24 13:08:05 +01:00
|
|
|
blas_int lda = ld;
|
2019-01-08 17:12:05 +01:00
|
|
|
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
|
2019-01-25 15:27:20 +01:00
|
|
|
ConstGeneralMatrix::multInvLeft(const std::string &trans, int mrows, int mcols, int mld, double *d) const
|
2019-01-08 17:12:05 +01:00
|
|
|
{
|
|
|
|
if (rows != cols)
|
2019-01-16 17:52:16 +01:00
|
|
|
throw SYLV_MES_EXCEPTION("The matrix is not square for inversion.");
|
2019-01-08 17:12:05 +01:00
|
|
|
if (cols != mrows)
|
2019-01-16 17:52:16 +01:00
|
|
|
throw SYLV_MES_EXCEPTION("Wrong dimensions for matrix inverse mutliply.");
|
2019-01-08 17:12:05 +01:00
|
|
|
|
|
|
|
if (rows > 0)
|
|
|
|
{
|
|
|
|
GeneralMatrix inv(*this);
|
2019-01-16 17:52:16 +01:00
|
|
|
std::vector<lapack_int> ipiv(rows);
|
2019-01-08 17:12:05 +01:00
|
|
|
lapack_int info;
|
2019-01-24 13:08:05 +01:00
|
|
|
lapack_int rows2 = rows, mcols2 = mcols, mld2 = mld, lda = inv.ld;
|
|
|
|
dgetrf(&rows2, &rows2, inv.getData().base(), &lda, ipiv.data(), &info);
|
2019-01-25 15:27:20 +01:00
|
|
|
dgetrs(trans.c_str(), &rows2, &mcols2, inv.base(), &lda, ipiv.data(), d,
|
2019-01-08 17:12:05 +01:00
|
|
|
&mld2, &info);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
/* 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)
|
2019-03-28 18:56:46 +01:00
|
|
|
throw SYLV_MES_EXCEPTION(u8"Skip≠1 not implemented in ConstGeneralMatrix::multInvLeft(Vector&)");
|
2019-01-08 17:12:05 +01:00
|
|
|
|
|
|
|
multInvLeft("N", d.length(), 1, d.length(), d.base());
|
|
|
|
}
|
|
|
|
|
|
|
|
/* d = inv(this')*d */
|
|
|
|
void
|
|
|
|
ConstGeneralMatrix::multInvLeftTrans(Vector &d) const
|
|
|
|
{
|
|
|
|
if (d.skip() != 1)
|
2019-03-28 18:56:46 +01:00
|
|
|
throw SYLV_MES_EXCEPTION(u8"Skip≠1 not implemented in ConstGeneralMatrix::multInvLeft(Vector&)");
|
2019-01-08 17:12:05 +01:00
|
|
|
|
|
|
|
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
|
|
|
|
{
|
2019-01-16 17:52:16 +01:00
|
|
|
auto ff = std::cout.flags();
|
|
|
|
std::cout << "rows=" << rows << ", cols=" << cols << std::endl;
|
2019-01-08 17:12:05 +01:00
|
|
|
for (int i = 0; i < rows; i++)
|
|
|
|
{
|
2019-01-16 17:52:16 +01:00
|
|
|
std::cout << "row " << i << ':' << std::endl
|
|
|
|
<< std::setprecision(3);
|
2019-01-08 17:12:05 +01:00
|
|
|
for (int j = 0; j < cols; j++)
|
2019-01-16 17:52:16 +01:00
|
|
|
std::cout << std::setw(6) << get(i, j) << ' ';
|
|
|
|
std::cout << std::endl;
|
2019-01-08 17:12:05 +01:00
|
|
|
}
|
2019-01-16 17:52:16 +01:00
|
|
|
std::cout.flags(ff);
|
2019-01-08 17:12:05 +01:00
|
|
|
}
|
|
|
|
|
|
|
|
void
|
|
|
|
SVDDecomp::construct(const GeneralMatrix &A)
|
|
|
|
{
|
|
|
|
// quick exit if empty matrix
|
|
|
|
if (minmn == 0)
|
|
|
|
{
|
|
|
|
U.unit();
|
|
|
|
VT.unit();
|
|
|
|
conv = true;
|
|
|
|
return;
|
|
|
|
}
|
|
|
|
|
|
|
|
// make copy of the matrix
|
|
|
|
GeneralMatrix AA(A);
|
|
|
|
|
|
|
|
lapack_int m = AA.numRows();
|
|
|
|
lapack_int n = AA.numCols();
|
|
|
|
double *a = AA.base();
|
|
|
|
lapack_int lda = AA.getLD();
|
|
|
|
double *s = sigma.base();
|
|
|
|
double *u = U.base();
|
|
|
|
lapack_int ldu = U.getLD();
|
|
|
|
double *vt = VT.base();
|
|
|
|
lapack_int ldvt = VT.getLD();
|
|
|
|
double tmpwork;
|
|
|
|
lapack_int lwork = -1;
|
|
|
|
lapack_int info;
|
|
|
|
|
2019-01-16 17:52:16 +01:00
|
|
|
std::vector<lapack_int> iwork(8*minmn);
|
2019-01-08 17:12:05 +01:00
|
|
|
// query for optimal lwork
|
|
|
|
dgesdd("A", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, &tmpwork,
|
2019-01-16 17:52:16 +01:00
|
|
|
&lwork, iwork.data(), &info);
|
2019-03-08 15:32:13 +01:00
|
|
|
lwork = static_cast<lapack_int>(tmpwork);
|
2019-01-08 17:12:05 +01:00
|
|
|
Vector work(lwork);
|
|
|
|
// do the decomposition
|
|
|
|
dgesdd("A", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work.base(),
|
2019-01-16 17:52:16 +01:00
|
|
|
&lwork, iwork.data(), &info);
|
2019-01-08 17:12:05 +01:00
|
|
|
if (info < 0)
|
|
|
|
throw SYLV_MES_EXCEPTION("Internal error in SVDDecomp constructor");
|
|
|
|
if (info == 0)
|
|
|
|
conv = true;
|
|
|
|
}
|
|
|
|
|
|
|
|
void
|
2019-01-24 15:22:36 +01:00
|
|
|
SVDDecomp::solve(const ConstGeneralMatrix &B, GeneralMatrix &X) const
|
2019-01-08 17:12:05 +01:00
|
|
|
{
|
|
|
|
if (B.numRows() != U.numRows())
|
|
|
|
throw SYLV_MES_EXCEPTION("Incompatible number of rows ");
|
|
|
|
|
2019-03-27 19:22:35 +01:00
|
|
|
/* Reciprocal condition number for determination of zeros in the
|
|
|
|
end of sigma */
|
2019-03-08 15:32:13 +01:00
|
|
|
constexpr double rcond = 1e-13;
|
2019-01-08 17:12:05 +01:00
|
|
|
|
|
|
|
// determine nz=number of zeros in the end of sigma
|
|
|
|
int nz = 0;
|
|
|
|
while (nz < minmn && sigma[minmn-1-nz] < rcond*sigma[0])
|
|
|
|
nz++;
|
|
|
|
// take relevant B for sigma inversion
|
|
|
|
int m = U.numRows();
|
|
|
|
int n = VT.numCols();
|
2019-02-26 16:54:34 +01:00
|
|
|
GeneralMatrix Bprime(transpose(U) * B, m-minmn, 0, minmn-nz, B.numCols());
|
2019-01-08 17:12:05 +01:00
|
|
|
// solve sigma
|
|
|
|
for (int i = 0; i < minmn-nz; i++)
|
Dynare++ / sylvester equation solver: refactor Vector and ConstVector classes
- these classes now encapsulate a std::shared_ptr<{const, }double>, so that
they do not perform memory management, and several {Const,}Vector instances
can transparently share the same underlying data
- make converting constructor from ConstVector to Vector explicit, since that
entails memory allocation (but the reverse conversion is almost costless, so
keep it implicit); do the same for GeneralMatrix/ConstGeneralMatrix,
TwoDMatrix/ConstTwoDMatrix
- remove the constructors that were extracting a row/column from a matrix, and
replace them by getRow() and getCol() methods on {Const,}GeneralMatrix
- rename and change the API of the complex version Vector::add(), so that it is
explicit that it deals with complex numbers
- add constructors that take a MATLAB mxArray
2019-01-22 16:07:44 +01:00
|
|
|
Bprime.getRow(i).mult(1.0/sigma[i]);
|
2019-01-08 17:12:05 +01:00
|
|
|
// solve VT
|
|
|
|
X.zeros();
|
|
|
|
//- copy Bprime to right place of X
|
|
|
|
for (int i = 0; i < minmn-nz; i++)
|
Dynare++ / sylvester equation solver: refactor Vector and ConstVector classes
- these classes now encapsulate a std::shared_ptr<{const, }double>, so that
they do not perform memory management, and several {Const,}Vector instances
can transparently share the same underlying data
- make converting constructor from ConstVector to Vector explicit, since that
entails memory allocation (but the reverse conversion is almost costless, so
keep it implicit); do the same for GeneralMatrix/ConstGeneralMatrix,
TwoDMatrix/ConstTwoDMatrix
- remove the constructors that were extracting a row/column from a matrix, and
replace them by getRow() and getCol() methods on {Const,}GeneralMatrix
- rename and change the API of the complex version Vector::add(), so that it is
explicit that it deals with complex numbers
- add constructors that take a MATLAB mxArray
2019-01-22 16:07:44 +01:00
|
|
|
X.getRow(n-minmn+i) = Bprime.getRow(i);
|
2019-01-08 17:12:05 +01:00
|
|
|
//- multiply with VT
|
|
|
|
X.multLeftTrans(VT);
|
|
|
|
}
|