Adding separate and modified sylv utilities needed by KF
git-svn-id: https://www.dynare.org/svn/dynare/trunk@2763 ac1d8469-bf42-47a9-8791-bf33cf982152time-shift
george
parent
f74c6f30b2
commit
6687ebb992
|
|
@ -0,0 +1,497 @@
|
|||
/* $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.h"
|
||||
#include "GeneralMatrix.h"
|
||||
|
||||
#include "cppblas.h"
|
||||
#include "cpplapack.h"
|
||||
|
||||
#include <stdio.h>
|
||||
#include <string.h>
|
||||
#include <stdlib.h>
|
||||
#include <cmath>
|
||||
#include <limits>
|
||||
|
||||
int GeneralMatrix::md_length = 32;
|
||||
|
||||
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, const char* dummy)
|
||||
: data(m.rows*m.cols), rows(m.cols), cols(m.rows), ld(m.cols)
|
||||
{
|
||||
for (int i = 0; i < m.rows; i++)
|
||||
for (int j = 0; j < m.cols; j++)
|
||||
get(j,i) = m.get(i,j);
|
||||
}
|
||||
|
||||
GeneralMatrix::GeneralMatrix(const ConstGeneralMatrix& m, const char* dummy)
|
||||
: data(m.rows*m.cols), rows(m.cols), cols(m.rows), ld(m.cols)
|
||||
{
|
||||
for (int i = 0; i < m.rows; i++)
|
||||
for (int j = 0; j < m.cols; j++)
|
||||
get(j,i) = m.get(i,j);
|
||||
}
|
||||
|
||||
|
||||
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)
|
||||
: data(m.base()+m.ld*j+i, m.ld*(ncols-1)+nrows), rows(nrows), cols(ncols), ld(m.ld)
|
||||
{}
|
||||
|
||||
GeneralMatrix::GeneralMatrix(const GeneralMatrix& a, const GeneralMatrix& b)
|
||||
: data(a.rows*b.cols), rows(a.rows), cols(b.cols), ld(a.rows)
|
||||
{
|
||||
gemm("N", a, "N", b, 1.0, 0.0);
|
||||
}
|
||||
|
||||
GeneralMatrix::GeneralMatrix(const GeneralMatrix& a, const GeneralMatrix& b, const char* dum)
|
||||
: data(a.rows*b.rows), rows(a.rows), cols(b.rows), ld(a.rows)
|
||||
{
|
||||
gemm("N", a, "T", b, 1.0, 0.0);
|
||||
}
|
||||
|
||||
GeneralMatrix::GeneralMatrix(const GeneralMatrix& a, const char* dum, const GeneralMatrix& b)
|
||||
: data(a.cols*b.cols), rows(a.cols), cols(b.cols), ld(a.cols)
|
||||
{
|
||||
gemm("T", a, "N", b, 1.0, 0.0);
|
||||
}
|
||||
|
||||
GeneralMatrix::GeneralMatrix(const GeneralMatrix& a, const char* dum1,
|
||||
const GeneralMatrix& b, const char* dum2)
|
||||
: data(a.cols*b.rows), rows(a.cols), cols(b.rows), ld(a.cols)
|
||||
{
|
||||
gemm("T", a, "T", b, 1.0, 0.0);
|
||||
}
|
||||
|
||||
|
||||
|
||||
GeneralMatrix::~GeneralMatrix()
|
||||
{
|
||||
}
|
||||
|
||||
|
||||
|
||||
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,
|
||||
const char* dum, double mult)
|
||||
{
|
||||
gemm("N", a, "T", b, mult, 1.0);
|
||||
}
|
||||
|
||||
void GeneralMatrix::multAndAdd(const ConstGeneralMatrix& a, const char* dum,
|
||||
const ConstGeneralMatrix& b, double mult)
|
||||
{
|
||||
gemm("T", a, "N", b, mult, 1.0);
|
||||
}
|
||||
|
||||
void GeneralMatrix::multAndAdd(const ConstGeneralMatrix& a, const char* dum1,
|
||||
const ConstGeneralMatrix& b, const char* dum2, double mult)
|
||||
{
|
||||
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 {
|
||||
for (int i = 0; i < rows; i++)
|
||||
for (int j = 0; j < cols; j++)
|
||||
get(i,j) = 0;
|
||||
}
|
||||
}
|
||||
|
||||
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 {
|
||||
for (int i = 0; i < rows; i++)
|
||||
for (int j = 0; j < cols; j++)
|
||||
get(i,j) *= a;
|
||||
}
|
||||
}
|
||||
|
||||
// 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 {
|
||||
for (int i = 0; i < rows; i++)
|
||||
for (int j = 0; j < cols; j++)
|
||||
get(i,j) += a*m.get(i,j);
|
||||
}
|
||||
}
|
||||
|
||||
void GeneralMatrix::add(double a, const ConstGeneralMatrix& m, const char* dum)
|
||||
{
|
||||
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);
|
||||
}
|
||||
|
||||
|
||||
inline void
|
||||
GeneralMatrix::copy(const ConstGeneralMatrix& m)
|
||||
{
|
||||
memcpy(data.base() ,m.getData().base() ,m.numCols()*m.numRows()*sizeof(double));
|
||||
};
|
||||
|
||||
|
||||
inline void
|
||||
GeneralMatrix::copy(const GeneralMatrix& m)
|
||||
{
|
||||
memcpy(data.base(),m.getData().base(),m.numCols()*m.numRows()*sizeof(double));
|
||||
};
|
||||
|
||||
|
||||
void GeneralMatrix::gemm(const char* transa, const ConstGeneralMatrix& a,
|
||||
const char* transb, const ConstGeneralMatrix& b,
|
||||
double alpha, double beta)
|
||||
{
|
||||
int opa_rows = a.numRows();
|
||||
int opa_cols = a.numCols();
|
||||
if (!strcmp(transa, "T")) {
|
||||
opa_rows = a.numCols();
|
||||
opa_cols = a.numRows();
|
||||
}
|
||||
int opb_rows = b.numRows();
|
||||
int opb_cols = b.numCols();
|
||||
if (!strcmp(transb, "T")) {
|
||||
opb_rows = b.numCols();
|
||||
opb_cols = b.numRows();
|
||||
}
|
||||
|
||||
if (opa_rows != numRows() ||
|
||||
opb_cols != numCols() ||
|
||||
opa_cols != opb_rows) {
|
||||
throw SYLV_MES_EXCEPTION("Wrong dimensions for matrix multiplication.");
|
||||
}
|
||||
|
||||
int m = opa_rows;
|
||||
int n = opb_cols;
|
||||
int k = opa_cols;
|
||||
int lda = a.ld;
|
||||
int ldb = b.ld;
|
||||
int ldc = ld;
|
||||
if (lda > 0 && ldb > 0 && ldc > 0) {
|
||||
BLAS_dgemm(transa, transb, &m, &n, &k, &alpha, a.data.base(), &lda,
|
||||
b.data.base(), &ldb, &beta, data.base(), &ldc);
|
||||
} else if (numRows()*numCols() > 0) {
|
||||
if (beta == 0.0)
|
||||
zeros();
|
||||
else
|
||||
mult(beta);
|
||||
}
|
||||
}
|
||||
|
||||
void GeneralMatrix::gemm_partial_left(const char* trans, const ConstGeneralMatrix& m,
|
||||
double alpha, double beta)
|
||||
{
|
||||
int icol;
|
||||
for (icol = 0; icol + md_length < cols; icol += md_length) {
|
||||
GeneralMatrix incopy((const GeneralMatrix&)*this, 0, icol, rows, md_length);
|
||||
GeneralMatrix inplace((GeneralMatrix&)*this, 0, icol, rows, md_length);
|
||||
inplace.gemm(trans, m, "N", ConstGeneralMatrix(incopy), alpha, beta);
|
||||
}
|
||||
if (cols > icol) {
|
||||
GeneralMatrix incopy((const GeneralMatrix&)*this, 0, icol, rows, cols - icol);
|
||||
GeneralMatrix inplace((GeneralMatrix&)*this, 0, icol, rows, cols - icol);
|
||||
inplace.gemm(trans, m, "N", ConstGeneralMatrix(incopy), alpha, beta);
|
||||
}
|
||||
}
|
||||
|
||||
void GeneralMatrix::gemm_partial_right(const char* trans, const ConstGeneralMatrix& m,
|
||||
double alpha, double beta)
|
||||
{
|
||||
int irow;
|
||||
for (irow = 0; irow + md_length < rows; irow += md_length) {
|
||||
GeneralMatrix incopy((const GeneralMatrix&)*this, irow, 0, md_length, cols);
|
||||
GeneralMatrix inplace((GeneralMatrix&)*this, irow, 0, md_length, cols);
|
||||
inplace.gemm("N", ConstGeneralMatrix(incopy), trans, m, alpha, beta);
|
||||
}
|
||||
if (rows > irow) {
|
||||
GeneralMatrix incopy((const GeneralMatrix&)*this, irow, 0, rows - irow, cols);
|
||||
GeneralMatrix inplace((GeneralMatrix&)*this, irow, 0, rows - irow, cols);
|
||||
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())
|
||||
{
|
||||
// can check that the submatirx is fully in the matrix
|
||||
}
|
||||
|
||||
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())
|
||||
{
|
||||
// can check that the submatirx is fully in the matrix
|
||||
}
|
||||
|
||||
|
||||
ConstGeneralMatrix::ConstGeneralMatrix(const GeneralMatrix& m)
|
||||
: data(m.data), rows(m.rows), cols(m.cols), ld(m.ld) {}
|
||||
|
||||
double ConstGeneralMatrix::getNormInf() const
|
||||
{
|
||||
double norm = 0.0;
|
||||
for (int i = 0; i < numRows(); i++) {
|
||||
ConstVector rowi(data.base()+i, ld, cols);
|
||||
double normi = rowi.getNorm1();
|
||||
if (norm < normi)
|
||||
norm = normi;
|
||||
}
|
||||
return norm;
|
||||
}
|
||||
|
||||
double ConstGeneralMatrix::getNorm1() const
|
||||
{
|
||||
double norm = 0.0;
|
||||
for (int j = 0; j < numCols(); j++) {
|
||||
ConstVector colj(data.base()+ld*j, 1, rows);
|
||||
double normj = colj.getNorm1();
|
||||
if (norm < normj)
|
||||
norm = normj;
|
||||
}
|
||||
return norm;
|
||||
}
|
||||
|
||||
void ConstGeneralMatrix::multVec(double a, Vector& x, double b, const ConstVector& d) const
|
||||
{
|
||||
if (x.length() != rows || cols != d.length()) {
|
||||
throw SYLV_MES_EXCEPTION("Wrong dimensions for vector multiply.");
|
||||
}
|
||||
if (rows > 0) {
|
||||
int mm = rows;
|
||||
int nn = cols;
|
||||
double alpha = b;
|
||||
int lda = ld;
|
||||
int incx = d.skip();
|
||||
double beta = a;
|
||||
int incy = x.skip();
|
||||
BLAS_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()) {
|
||||
throw SYLV_MES_EXCEPTION("Wrong dimensions for vector multiply.");
|
||||
}
|
||||
if (rows > 0) {
|
||||
int mm = rows;
|
||||
int nn = cols;
|
||||
double alpha = b;
|
||||
int lda = rows;
|
||||
int incx = d.skip();
|
||||
double beta = a;
|
||||
int incy = x.skip();
|
||||
BLAS_dgemv("T", &mm, &nn, &alpha, data.base(), &lda, d.base(), &incx,
|
||||
&beta, x.base(), &incy);
|
||||
}
|
||||
}
|
||||
|
||||
/* m = inv(this)*m */
|
||||
void ConstGeneralMatrix::multInvLeft(const char* trans, int mrows, int mcols, int mld, double* d) const
|
||||
{
|
||||
if (rows != cols) {
|
||||
throw SYLV_MES_EXCEPTION("The matrix is not square for inversion.");
|
||||
}
|
||||
if (cols != mrows) {
|
||||
throw SYLV_MES_EXCEPTION("Wrong dimensions for matrix inverse mutliply.");
|
||||
}
|
||||
|
||||
if (rows > 0) {
|
||||
GeneralMatrix inv(*this);
|
||||
int* ipiv = new int[rows];
|
||||
int info;
|
||||
LAPACK_dgetrf(&rows, &rows, inv.getData().base(), &rows, ipiv, &info);
|
||||
LAPACK_dgetrs(trans, &rows, &mcols, inv.base(), &rows, ipiv, d,
|
||||
&mld, &info);
|
||||
delete [] ipiv;
|
||||
}
|
||||
}
|
||||
|
||||
/* 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) {
|
||||
throw SYLV_MES_EXCEPTION("Skip!=1 not implemented in ConstGeneralMatrix::multInvLeft(Vector&)");
|
||||
}
|
||||
|
||||
multInvLeft("N", d.length(), 1, d.length(), d.base());
|
||||
}
|
||||
|
||||
/* d = inv(this')*d */
|
||||
void ConstGeneralMatrix::multInvLeftTrans(Vector& d) const
|
||||
{
|
||||
if (d.skip() != 1) {
|
||||
throw SYLV_MES_EXCEPTION("Skip!=1 not implemented in ConstGeneralMatrix::multInvLeft(Vector&)");
|
||||
}
|
||||
|
||||
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
|
||||
{
|
||||
printf("rows=%d, cols=%d\n",rows, cols);
|
||||
for (int i = 0; i < rows; i++) {
|
||||
printf("row %d:\n",i);
|
||||
for (int j = 0; j < cols; j++) {
|
||||
printf("%6.3g ",get(i,j));
|
||||
}
|
||||
printf("\n");
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,289 @@
|
|||
/* $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)
|
||||
void copy(const ConstGeneralMatrix& m, int ioff , int joff );
|
||||
void copy(const GeneralMatrix& m, int ioff , int joff )
|
||||
{copy(ConstGeneralMatrix(m), ioff, joff);}
|
||||
|
||||
void copy(const ConstGeneralMatrix& m);
|
||||
void copy(const GeneralMatrix& m);
|
||||
|
||||
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:
|
||||
|
|
@ -0,0 +1,41 @@
|
|||
# $Header: /var/lib/cvs/dynare_cpp/sylv/cc/Makefile,v 1.4 2005/01/18 21:28:26 kamenik Exp $
|
||||
|
||||
# Tag $Name: $
|
||||
|
||||
CC_FLAGS := -Wall -I../testing -I../cc -DMATLAB \
|
||||
$(CC_INCLUDE_PATH) -Ic:/"Program Files"/MATLAB_SV71/extern/include
|
||||
|
||||
LDFLAGS = -Wl,-L"c:/Program Files"/MATLAB_SV71/extern/lib/win32/microsoft/ \
|
||||
-Wl,-llibmex -Wl,-llibmx -Wl,-llibmwlapack -Wl,-llibdflapack \
|
||||
-lg2c -lmingw32 -lstdc++
|
||||
|
||||
LD_LIBS=$(LDFLAGS)
|
||||
|
||||
ifeq ($(DEBUG),yes)
|
||||
# CC_FLAGS := $(CC_FLAGS) -g -DTL_DEBUG=2
|
||||
CC_FLAGS := $(CC_FLAGS) -g -DPOSIX_THREADS
|
||||
else
|
||||
CC_FLAGS := $(CC_FLAGS) -O2
|
||||
endif
|
||||
|
||||
ifeq ($(OS),Windows_NT)
|
||||
CC_FLAGS := -mno-cygwin -mthreads $(CC_FLAGS)
|
||||
endif
|
||||
|
||||
|
||||
objects := $(patsubst %.cpp,%.o,$(wildcard *.cpp))
|
||||
headers := $(wildcard *.h)
|
||||
|
||||
all: $(objects)
|
||||
|
||||
sylvester.a: $(objects)
|
||||
ar cr sylvester.a $(objects)
|
||||
ranlib sylvester.a
|
||||
|
||||
clear:
|
||||
rm -f *.o
|
||||
rm -f sylvester.a
|
||||
|
||||
%.o : %.cpp $(headers)
|
||||
$(CC) $(CC_FLAGS) $(EXTERN_DEFS) -c $*.cpp
|
||||
|
||||
|
|
@ -0,0 +1,39 @@
|
|||
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/SylvException.h,v 1.1.1.1 2004/06/04 13:00:44 kamenik Exp $ */
|
||||
|
||||
/* Tag $Name: $ */
|
||||
|
||||
#ifndef SYLV_EXCEPTION_H
|
||||
#define SYLV_EXCEPTION_H
|
||||
|
||||
#include "SylvMemory.h"
|
||||
|
||||
|
||||
class SylvException : public MallocAllocator {
|
||||
protected:
|
||||
char file[50];
|
||||
int line;
|
||||
const SylvException* source;
|
||||
public:
|
||||
SylvException(const char* f, int l, const SylvException* s);
|
||||
virtual ~SylvException();
|
||||
virtual int printMessage(char* str, int maxlen) const;
|
||||
void printMessage() const;
|
||||
};
|
||||
|
||||
class SylvExceptionMessage : public SylvException {
|
||||
char message[500];
|
||||
public:
|
||||
SylvExceptionMessage(const char* f, int l, const char* mes);
|
||||
virtual int printMessage(char* str, int maxlen) const;
|
||||
};
|
||||
|
||||
// define macros:
|
||||
#define SYLV_EXCEPTION(exc) (SylvException(__FILE__, __LINE__, exc))
|
||||
#define SYLV_MES_EXCEPTION(mes) (SylvExceptionMessage(__FILE__, __LINE__, mes))
|
||||
|
||||
#endif /* SYLV_EXCEPTION_H */
|
||||
|
||||
|
||||
// Local Variables:
|
||||
// mode:C++
|
||||
// End:
|
||||
|
|
@ -0,0 +1,366 @@
|
|||
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/Vector.cpp,v 1.1.1.1 2004/06/04 13:01:13 kamenik Exp $ */
|
||||
|
||||
/* Tag $Name: $ */
|
||||
|
||||
|
||||
#include "Vector.h"
|
||||
#include "GeneralMatrix.h"
|
||||
#include "SylvException.h"
|
||||
#include "cppblas.h"
|
||||
|
||||
#include <stdio.h>
|
||||
#include <string.h>
|
||||
#include <stdlib.h>
|
||||
#include <cmath>
|
||||
#include <algorithm>
|
||||
#include <limits>
|
||||
|
||||
using namespace std;
|
||||
|
||||
ZeroPad zero_pad;
|
||||
|
||||
Vector::Vector(const ConstVector& v)
|
||||
: len(v.length()), s(1), data(new double[len]), destroy(true)
|
||||
{
|
||||
copy(v.base(), v.skip());
|
||||
}
|
||||
|
||||
const Vector& Vector::operator=(const Vector& v)
|
||||
{
|
||||
if (this == &v)
|
||||
return *this;
|
||||
|
||||
if (v.length() != length()) {
|
||||
throw SYLV_MES_EXCEPTION("Attempt to assign vectors with different lengths.");
|
||||
}
|
||||
if (s == v.s &&
|
||||
(data <= v.data && v.data < data+len*s ||
|
||||
v.data <= data && data < v.data+v.len*v.s) &&
|
||||
(data-v.data) % s == 0) {
|
||||
printf("this destroy=%d, v destroy=%d, data-v.data=%d, len=%d\n", destroy, v.destroy, data-v.data, len);
|
||||
throw SYLV_MES_EXCEPTION("Attempt to assign overlapping vectors.");
|
||||
}
|
||||
copy(v.base(), v.skip());
|
||||
return *this;
|
||||
}
|
||||
|
||||
const Vector& Vector::operator=(const ConstVector& v)
|
||||
{
|
||||
if (v.length() != length()) {
|
||||
throw SYLV_MES_EXCEPTION("Attempt to assign vectors with different lengths.");
|
||||
}
|
||||
if (v.skip() == 1 && skip() == 1 && (
|
||||
(base() < v.base() + v.length() && base() >= v.base()) ||
|
||||
(base() + length() < v.base() + v.length() &&
|
||||
base() + length() > v.base()))) {
|
||||
throw SYLV_MES_EXCEPTION("Attempt to assign overlapping vectors.");
|
||||
}
|
||||
copy(v.base(), v.skip());
|
||||
return *this;
|
||||
}
|
||||
|
||||
void Vector::copy(const double* d, int inc)
|
||||
{
|
||||
int n = length();
|
||||
int incy = skip();
|
||||
BLAS_dcopy(&n, d, &inc, base(), &incy);
|
||||
}
|
||||
|
||||
Vector::Vector(Vector& v, int off, int l)
|
||||
: len(l), s(v.skip()), data(v.base()+off*v.skip()), destroy(false)
|
||||
{
|
||||
if (off < 0 || off + length() > v.length())
|
||||
throw SYLV_MES_EXCEPTION("Subvector not contained in supvector.");
|
||||
}
|
||||
|
||||
Vector::Vector(const Vector& v, int off, int l)
|
||||
: len(l), s(1), data(new double[len]), destroy(true)
|
||||
{
|
||||
if (off < 0 || off + length() > v.length())
|
||||
throw SYLV_MES_EXCEPTION("Subvector not contained in supvector.");
|
||||
copy(v.base()+off*v.skip(), v.skip());
|
||||
}
|
||||
|
||||
Vector::Vector(GeneralMatrix& m, int col)
|
||||
: len(m.numRows()), s(1), data(&(m.get(0, col))), destroy(false)
|
||||
{
|
||||
}
|
||||
|
||||
Vector::Vector(int row, GeneralMatrix& m)
|
||||
: len(m.numCols()), s(m.getLD()), data(&(m.get(row, 0))), destroy(false)
|
||||
{
|
||||
}
|
||||
|
||||
bool Vector::operator==(const Vector& y) const
|
||||
{
|
||||
return ConstVector(*this) == y;
|
||||
}
|
||||
|
||||
bool Vector::operator!=(const Vector& y) const
|
||||
{
|
||||
return ConstVector(*this) != y;
|
||||
}
|
||||
|
||||
bool Vector::operator<(const Vector& y) const
|
||||
{
|
||||
return ConstVector(*this) < y;
|
||||
}
|
||||
|
||||
bool Vector::operator<=(const Vector& y) const
|
||||
{
|
||||
return ConstVector(*this) <= y;
|
||||
}
|
||||
|
||||
bool Vector::operator>(const Vector& y) const
|
||||
{
|
||||
return ConstVector(*this) > y;
|
||||
}
|
||||
|
||||
bool Vector::operator>=(const Vector& y) const
|
||||
{
|
||||
return ConstVector(*this) >= y;
|
||||
}
|
||||
|
||||
void Vector::zeros()
|
||||
{
|
||||
if (skip() == 1) {
|
||||
double* p = base();
|
||||
for (int i = 0; i < length()/ZeroPad::length;
|
||||
i++, p += ZeroPad::length)
|
||||
memcpy(p, zero_pad.getBase(), sizeof(double)*ZeroPad::length);
|
||||
for ( ; p < base()+length(); p++)
|
||||
*p = 0.0;
|
||||
} else {
|
||||
for (int i = 0; i < length(); i++)
|
||||
operator[](i) = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
void Vector::nans()
|
||||
{
|
||||
for (int i = 0; i < length(); i++)
|
||||
operator[](i) = std::numeric_limits<double>::quiet_NaN();
|
||||
}
|
||||
|
||||
void Vector::infs()
|
||||
{
|
||||
for (int i = 0; i < length(); i++)
|
||||
operator[](i) = std::numeric_limits<double>::infinity();
|
||||
}
|
||||
|
||||
Vector::~Vector()
|
||||
{
|
||||
if (destroy) {
|
||||
delete [] data;
|
||||
}
|
||||
}
|
||||
|
||||
void Vector::rotatePair(double alpha, double beta1, double beta2, int i)
|
||||
{
|
||||
double tmp = alpha*operator[](i) - beta1*operator[](i+1);
|
||||
operator[](i+1) = alpha*operator[](i+1) - beta2*operator[](i);
|
||||
operator[](i) = tmp;
|
||||
}
|
||||
|
||||
void Vector::add(double r, const Vector& v)
|
||||
{
|
||||
add(r, ConstVector(v));
|
||||
}
|
||||
|
||||
void Vector::add(double r, const ConstVector& v)
|
||||
{
|
||||
int n = length();
|
||||
int incx = v.skip();
|
||||
int incy = skip();
|
||||
BLAS_daxpy(&n, &r, v.base(), &incx, base(), &incy);
|
||||
}
|
||||
|
||||
void Vector::add(const double* z, const Vector& v)
|
||||
{
|
||||
add(z, ConstVector(v));
|
||||
}
|
||||
|
||||
void Vector::add(const double* z, const ConstVector& v)
|
||||
{
|
||||
int n = length()/2;
|
||||
int incx = v.skip();
|
||||
int incy = skip();
|
||||
BLAS_zaxpy(&n, z, v.base(), &incx, base(), &incy);
|
||||
}
|
||||
|
||||
void Vector::mult(double r)
|
||||
{
|
||||
int n = length();
|
||||
int incx = skip();
|
||||
BLAS_dscal(&n, &r, base(), &incx);
|
||||
}
|
||||
|
||||
void Vector::mult2(double alpha, double beta1, double beta2,
|
||||
Vector& x1, Vector& x2,
|
||||
const Vector& b1, const Vector& b2)
|
||||
{
|
||||
x1.zeros();
|
||||
x2.zeros();
|
||||
mult2a(alpha, beta1, beta2, x1, x2, b1, b2);
|
||||
}
|
||||
|
||||
void Vector::mult2a(double alpha, double beta1, double beta2,
|
||||
Vector& x1, Vector& x2,
|
||||
const Vector& b1, const Vector& b2)
|
||||
{
|
||||
x1.add(alpha, b1);
|
||||
x1.add(-beta1, b2);
|
||||
x2.add(alpha, b2);
|
||||
x2.add(-beta2, b1);
|
||||
}
|
||||
|
||||
double Vector::getNorm() const
|
||||
{
|
||||
ConstVector v(*this);
|
||||
return v.getNorm();
|
||||
}
|
||||
|
||||
double Vector::getMax() const
|
||||
{
|
||||
ConstVector v(*this);
|
||||
return v.getMax();
|
||||
}
|
||||
|
||||
double Vector::getNorm1() const
|
||||
{
|
||||
ConstVector v(*this);
|
||||
return v.getNorm1();
|
||||
}
|
||||
|
||||
double Vector::dot(const Vector& y) const
|
||||
{
|
||||
return ConstVector(*this).dot(ConstVector(y));
|
||||
}
|
||||
|
||||
bool Vector::isFinite() const
|
||||
{
|
||||
return (ConstVector(*this)).isFinite();
|
||||
}
|
||||
|
||||
void Vector::print() const
|
||||
{
|
||||
for (int i = 0; i < length(); i++) {
|
||||
printf("%d\t%8.4g\n", i, operator[](i));
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
ConstVector::ConstVector(const Vector& v, int off, int l)
|
||||
: BaseConstVector(l, v.skip(), v.base() + v.skip()*off)
|
||||
{
|
||||
if (off < 0 || off + length() > v.length()) {
|
||||
throw SYLV_MES_EXCEPTION("Subvector not contained in supvector.");
|
||||
}
|
||||
}
|
||||
|
||||
ConstVector::ConstVector(const ConstVector& v, int off, int l)
|
||||
: BaseConstVector(l, v.skip(), v.base() + v.skip()*off)
|
||||
{
|
||||
if (off < 0 || off + length() > v.length()) {
|
||||
throw SYLV_MES_EXCEPTION("Subvector not contained in supvector.");
|
||||
}
|
||||
}
|
||||
|
||||
ConstVector::ConstVector(const double* d, int skip, int l)
|
||||
: BaseConstVector(l, skip, d)
|
||||
{
|
||||
}
|
||||
|
||||
ConstVector::ConstVector(const ConstGeneralMatrix& m, int col)
|
||||
: BaseConstVector(m.numRows(), 1, &(m.get(0, col)))
|
||||
{
|
||||
}
|
||||
|
||||
ConstVector::ConstVector(int row, const ConstGeneralMatrix& m)
|
||||
: BaseConstVector(m.numCols(), m.getLD(), &(m.get(row, 0)))
|
||||
{
|
||||
}
|
||||
|
||||
bool ConstVector::operator==(const ConstVector& y) const
|
||||
{
|
||||
if (length() != y.length())
|
||||
return false;
|
||||
if (length() == 0)
|
||||
return true;
|
||||
int i = 0;
|
||||
while (i < length() && operator[](i) == y[i])
|
||||
i++;
|
||||
return i == length();
|
||||
}
|
||||
|
||||
bool ConstVector::operator<(const ConstVector& y) const
|
||||
{
|
||||
int i = std::min(length(), y.length());
|
||||
int ii = 0;
|
||||
while (ii < i && operator[](ii) == y[ii])
|
||||
ii++;
|
||||
if (ii < i)
|
||||
return operator[](ii) < y[ii];
|
||||
else
|
||||
return length() < y.length();
|
||||
}
|
||||
|
||||
double ConstVector::getNorm() const
|
||||
{
|
||||
double s = 0;
|
||||
for (int i = 0; i < length(); i++) {
|
||||
s+=operator[](i)*operator[](i);
|
||||
}
|
||||
return sqrt(s);
|
||||
}
|
||||
|
||||
double ConstVector::getMax() const
|
||||
{
|
||||
double r = 0;
|
||||
for (int i = 0; i < length(); i++) {
|
||||
if (abs(operator[](i))>r)
|
||||
r = abs(operator[](i));
|
||||
}
|
||||
return r;
|
||||
}
|
||||
|
||||
double ConstVector::getNorm1() const
|
||||
{
|
||||
double norm = 0.0;
|
||||
for (int i = 0; i < length(); i++) {
|
||||
norm += abs(operator[](i));
|
||||
}
|
||||
return norm;
|
||||
}
|
||||
|
||||
double ConstVector::dot(const ConstVector& y) const
|
||||
{
|
||||
if (length() != y.length())
|
||||
throw SYLV_MES_EXCEPTION("Vector has different length in ConstVector::dot.");
|
||||
int n = length();
|
||||
int incx = skip();
|
||||
int incy = y.skip();
|
||||
return BLAS_ddot(&n, base(), &incx, y.base(), &incy);
|
||||
}
|
||||
|
||||
bool ConstVector::isFinite() const
|
||||
{
|
||||
int i = 0;
|
||||
while (i < length() && isfinite(operator[](i)))
|
||||
i++;
|
||||
return i == length();
|
||||
}
|
||||
|
||||
void ConstVector::print() const
|
||||
{
|
||||
for (int i = 0; i < length(); i++) {
|
||||
printf("%d\t%8.4g\n", i, operator[](i));
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
ZeroPad::ZeroPad()
|
||||
{
|
||||
for (int i = 0; i < length; i++)
|
||||
pad[i] = 0.0;
|
||||
}
|
||||
|
|
@ -0,0 +1,178 @@
|
|||
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/Vector.h,v 1.1.1.1 2004/06/04 13:01:13 kamenik Exp $ */
|
||||
|
||||
/* Tag $Name: $ */
|
||||
|
||||
#ifndef VECTOR_H
|
||||
#define VECTOR_H
|
||||
|
||||
/* NOTE! Vector and ConstVector have not common super class in order
|
||||
* to avoid running virtual method invokation mechanism. Some
|
||||
* members, and methods are thus duplicated */
|
||||
|
||||
#include <stdio.h>
|
||||
|
||||
class GeneralMatrix;
|
||||
class ConstVector;
|
||||
|
||||
class Vector {
|
||||
protected:
|
||||
int len;
|
||||
int s;
|
||||
double* data;
|
||||
bool destroy;
|
||||
public:
|
||||
Vector() : len(0), s(1), data(0), destroy(false) {}
|
||||
Vector(int l) : len(l), s(1), data(new double[l]), destroy(true) {}
|
||||
Vector(Vector& v) : len(v.length()), s(v.skip()), data(v.base()), destroy(false) {}
|
||||
Vector(const Vector& v)
|
||||
: len(v.length()), s(1), data(new double[len]), destroy(true)
|
||||
{copy(v.base(), v.skip());}
|
||||
Vector(const ConstVector& v);
|
||||
Vector(const double* d, int l)
|
||||
: len(l), s(1), data(new double[len]), destroy(true)
|
||||
{copy(d, 1);}
|
||||
Vector(double* d, int l)
|
||||
: len(l), s(1), data(d), destroy(false) {}
|
||||
Vector(double* d, int skip, int l)
|
||||
: len(l), s(skip), data(d), destroy(false) {}
|
||||
Vector(Vector& v, int off, int l);
|
||||
Vector(const Vector& v, int off, int l);
|
||||
Vector(GeneralMatrix& m, int col);
|
||||
Vector(int row, GeneralMatrix& m);
|
||||
const Vector& operator=(const Vector& v);
|
||||
const Vector& operator=(const ConstVector& v);
|
||||
double& operator[](int i)
|
||||
{return data[s*i];}
|
||||
const double& operator[](int i) const
|
||||
{return data[s*i];}
|
||||
const double* base() const
|
||||
{return data;}
|
||||
double* base()
|
||||
{return data;}
|
||||
int length() const
|
||||
{return len;}
|
||||
int skip() const
|
||||
{return s;}
|
||||
|
||||
/** Exact equality. */
|
||||
bool operator==(const Vector& y) const;
|
||||
bool operator!=(const Vector& y) const;
|
||||
/** Lexicographic ordering. */
|
||||
bool operator<(const Vector& y) const;
|
||||
bool operator<=(const Vector& y) const;
|
||||
bool operator>(const Vector& y) const;
|
||||
bool operator>=(const Vector& y) const;
|
||||
|
||||
virtual ~Vector();
|
||||
void zeros();
|
||||
void nans();
|
||||
void infs();
|
||||
bool toBeDestroyed() const {return destroy;}
|
||||
void rotatePair(double alpha, double beta1, double beta2, int i);
|
||||
void add(double r, const Vector& v);
|
||||
void add(double r, const ConstVector& v);
|
||||
void add(const double* z, const Vector& v);
|
||||
void add(const double* z, const ConstVector& v);
|
||||
void mult(double r);
|
||||
double getNorm() const;
|
||||
double getMax() const;
|
||||
double getNorm1() const;
|
||||
double dot(const Vector& y) const;
|
||||
bool isFinite() const;
|
||||
void print() const;
|
||||
|
||||
/* multiplies | alpha -beta1| |b1| |x1|
|
||||
| |\otimes I .| | = | |
|
||||
| -beta2 alpha| |b2| |x2|
|
||||
*/
|
||||
static void mult2(double alpha, double beta1, double beta2,
|
||||
Vector& x1, Vector& x2,
|
||||
const Vector& b1, const Vector& b2);
|
||||
/* same, but adds instead of set */
|
||||
static void mult2a(double alpha, double beta1, double beta2,
|
||||
Vector& x1, Vector& x2,
|
||||
const Vector& b1, const Vector& b2);
|
||||
/* same, but subtracts instead of add */
|
||||
static void mult2s(double alpha, double beta1, double beta2,
|
||||
Vector& x1, Vector& x2,
|
||||
const Vector& b1, const Vector& b2)
|
||||
{mult2a(-alpha, -beta1, -beta2, x1, x2, b1, b2);}
|
||||
private:
|
||||
void copy(const double* d, int inc);
|
||||
const Vector& operator=(int); // must not be used (not implemented)
|
||||
const Vector& operator=(double); // must not be used (not implemented)
|
||||
};
|
||||
|
||||
|
||||
class BaseConstVector {
|
||||
protected:
|
||||
int len;
|
||||
int s;
|
||||
const double* data;
|
||||
public:
|
||||
BaseConstVector(int l, int si, const double* d) : len(l), s(si), data(d) {}
|
||||
BaseConstVector(const BaseConstVector& v) : len(v.len), s(v.s), data(v.data) {}
|
||||
const BaseConstVector& operator=(const BaseConstVector& v)
|
||||
{len = v.len; s = v.s; data = v.data; return *this;}
|
||||
const double& operator[](int i) const
|
||||
{return data[s*i];}
|
||||
const double* base() const
|
||||
{return data;}
|
||||
int length() const
|
||||
{return len;}
|
||||
int skip() const
|
||||
{return s;}
|
||||
};
|
||||
|
||||
class ConstGeneralMatrix;
|
||||
|
||||
class ConstVector : public BaseConstVector {
|
||||
public:
|
||||
ConstVector(const Vector& v) : BaseConstVector(v.length(), v.skip(), v.base()) {}
|
||||
ConstVector(const ConstVector& v) : BaseConstVector(v) {}
|
||||
ConstVector(const double* d, int l) : BaseConstVector(l, 1, d) {}
|
||||
ConstVector(const Vector& v, int off, int l);
|
||||
ConstVector(const ConstVector& v, int off, int l);
|
||||
ConstVector(const double* d, int skip, int l);
|
||||
ConstVector(const ConstGeneralMatrix& m, int col);
|
||||
ConstVector(int row, const ConstGeneralMatrix& m);
|
||||
|
||||
virtual ~ConstVector() {}
|
||||
/** Exact equality. */
|
||||
bool operator==(const ConstVector& y) const;
|
||||
bool operator!=(const ConstVector& y) const
|
||||
{return ! operator==(y);}
|
||||
/** Lexicographic ordering. */
|
||||
bool operator<(const ConstVector& y) const;
|
||||
bool operator<=(const ConstVector& y) const
|
||||
{return operator<(y) || operator==(y);}
|
||||
bool operator>(const ConstVector& y) const
|
||||
{return ! operator<=(y);}
|
||||
bool operator>=(const ConstVector& y) const
|
||||
{return ! operator<(y);}
|
||||
|
||||
double getNorm() const;
|
||||
double getMax() const;
|
||||
double getNorm1() const;
|
||||
double dot(const ConstVector& y) const;
|
||||
bool isFinite() const;
|
||||
void print() const;
|
||||
};
|
||||
|
||||
class ZeroPad {
|
||||
public:
|
||||
//static const int length = 16;
|
||||
enum { length = 16};
|
||||
private:
|
||||
double pad[16];
|
||||
public:
|
||||
ZeroPad();
|
||||
const double* getBase() const {return pad;}
|
||||
};
|
||||
|
||||
#endif /* VECTOR_H */
|
||||
|
||||
|
||||
// Local Variables:
|
||||
// mode:C++
|
||||
// End:
|
||||
|
|
@ -0,0 +1,68 @@
|
|||
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/cppblas.h,v 1.2 2004/11/24 20:42:52 kamenik Exp $ */
|
||||
|
||||
/* Tag $Name: $ */
|
||||
|
||||
#ifndef CPPBLAS_H
|
||||
#define CPPBLAS_H
|
||||
|
||||
#if defined(MATLAB) && !defined(__linux__)
|
||||
#define BLAS_dgemm dgemm
|
||||
#define BLAS_dgemv dgemv
|
||||
#define BLAS_dtrsv dtrsv
|
||||
#define BLAS_dtrmv dtrmv
|
||||
#define BLAS_daxpy daxpy
|
||||
#define BLAS_dcopy dcopy
|
||||
#define BLAS_zaxpy zaxpy
|
||||
#define BLAS_dscal dscal
|
||||
#define BLAS_dtrsm dtrsm
|
||||
#define BLAS_ddot ddot
|
||||
#else /* defined(MATLAB) && !defined(__linux__) */
|
||||
#define BLAS_dgemm dgemm_
|
||||
#define BLAS_dgemv dgemv_
|
||||
#define BLAS_dtrsv dtrsv_
|
||||
#define BLAS_dtrmv dtrmv_
|
||||
#define BLAS_daxpy daxpy_
|
||||
#define BLAS_dcopy dcopy_
|
||||
#define BLAS_zaxpy zaxpy_
|
||||
#define BLAS_dscal dscal_
|
||||
#define BLAS_dtrsm dtrsm_
|
||||
#define BLAS_ddot ddot_
|
||||
#endif /* defined(MATLAB) && !defined(__linux__) */
|
||||
|
||||
#define BLCHAR const char*
|
||||
#define CONST_BLINT const int*
|
||||
#define CONST_BLDOU const double*
|
||||
#define BLDOU double*
|
||||
|
||||
extern "C" {
|
||||
void BLAS_dgemm(BLCHAR transa, BLCHAR transb, CONST_BLINT m, CONST_BLINT n,
|
||||
CONST_BLINT k, CONST_BLDOU alpha, CONST_BLDOU a, CONST_BLINT lda,
|
||||
CONST_BLDOU b, CONST_BLINT ldb, CONST_BLDOU beta,
|
||||
BLDOU c, CONST_BLINT ldc);
|
||||
void BLAS_dgemv(BLCHAR trans, CONST_BLINT m, CONST_BLINT n, CONST_BLDOU alpha,
|
||||
CONST_BLDOU a, CONST_BLINT lda, CONST_BLDOU x, CONST_BLINT incx,
|
||||
CONST_BLDOU beta, BLDOU y, CONST_BLINT incy);
|
||||
void BLAS_dtrsv(BLCHAR uplo, BLCHAR trans, BLCHAR diag, CONST_BLINT n,
|
||||
CONST_BLDOU a, CONST_BLINT lda, BLDOU x, CONST_BLINT incx);
|
||||
void BLAS_dtrmv(BLCHAR uplo, BLCHAR trans, BLCHAR diag, CONST_BLINT n,
|
||||
CONST_BLDOU a, CONST_BLINT lda, BLDOU x, CONST_BLINT incx);
|
||||
void BLAS_daxpy(CONST_BLINT n, CONST_BLDOU a, CONST_BLDOU x, CONST_BLINT incx,
|
||||
BLDOU y, CONST_BLINT incy);
|
||||
void BLAS_dcopy(CONST_BLINT n, CONST_BLDOU x, CONST_BLINT incx,
|
||||
BLDOU y, CONST_BLINT incy);
|
||||
void BLAS_zaxpy(CONST_BLINT n, CONST_BLDOU a, CONST_BLDOU x, CONST_BLINT incx,
|
||||
BLDOU y, CONST_BLINT incy);
|
||||
void BLAS_dscal(CONST_BLINT n, CONST_BLDOU a, BLDOU x, CONST_BLINT incx);
|
||||
void BLAS_dtrsm(BLCHAR side, BLCHAR uplo, BLCHAR transa, BLCHAR diag, CONST_BLINT m,
|
||||
CONST_BLINT n, CONST_BLDOU alpha, CONST_BLDOU a, CONST_BLINT lda,
|
||||
BLDOU b, CONST_BLINT ldb);
|
||||
double BLAS_ddot(CONST_BLINT n, CONST_BLDOU x, CONST_BLINT incx, CONST_BLDOU y,
|
||||
CONST_BLINT incy);
|
||||
};
|
||||
|
||||
|
||||
#endif /* CPPBLAS_H */
|
||||
|
||||
// Local Variables:
|
||||
// mode:C++
|
||||
// End:
|
||||
|
|
@ -0,0 +1,82 @@
|
|||
/* $Header: /var/lib/cvs/dynare_cpp/sylv/cc/cpplapack.h,v 1.3 2004/11/24 20:43:10 kamenik Exp $ */
|
||||
|
||||
/* Tag $Name: $ */
|
||||
|
||||
#ifndef CPPLAPACK_H
|
||||
#define CPPLAPACK_H
|
||||
|
||||
#if defined(MATLAB) && !defined(__linux__)
|
||||
#define LAPACK_dgetrs dgetrs
|
||||
#define LAPACK_dgetrf dgetrf
|
||||
#define LAPACK_dgees dgees
|
||||
#define LAPACK_dgecon dgecon
|
||||
#define LAPACK_dtrexc dtrexc
|
||||
#define LAPACK_dtrsyl dtrsyl
|
||||
#define LAPACK_dpotrf dpotrf
|
||||
#define LAPACK_dgges dgges
|
||||
#define LAPACK_dsyev dsyev
|
||||
#define LAPACK_dsyevr dsyevr
|
||||
#else /* MATLAB */
|
||||
#define LAPACK_dgetrs dgetrs_
|
||||
#define LAPACK_dgetrf dgetrf_
|
||||
#define LAPACK_dgees dgees_
|
||||
#define LAPACK_dgecon dgecon_
|
||||
#define LAPACK_dtrexc dtrexc_
|
||||
#define LAPACK_dtrsyl dtrsyl_
|
||||
#define LAPACK_dpotrf dpotrf_
|
||||
#define LAPACK_dgges dgges_
|
||||
#define LAPACK_dsyev dsyev_
|
||||
#define LAPACK_dsyevr dsyevr_
|
||||
#endif /* MATLAB */
|
||||
|
||||
#define LACHAR const char*
|
||||
#define CONST_LAINT const int*
|
||||
#define LAINT int*
|
||||
#define CONST_LADOU const double*
|
||||
#define LADOU double*
|
||||
typedef int (*DGGESCRIT)(const double*, const double*, const double*);
|
||||
|
||||
extern "C" {
|
||||
void LAPACK_dgetrs(LACHAR trans, CONST_LAINT n, CONST_LAINT nrhs,
|
||||
CONST_LADOU a, CONST_LAINT lda, CONST_LAINT ipiv,
|
||||
LADOU b, CONST_LAINT ldb, LAINT info);
|
||||
void LAPACK_dgetrf(CONST_LAINT m, CONST_LAINT n, LADOU a,
|
||||
CONST_LAINT lda, LAINT ipiv, LAINT info);
|
||||
void LAPACK_dgees(LACHAR jobvs, LACHAR sort, const void* select,
|
||||
CONST_LAINT n, LADOU a, CONST_LAINT lda, LAINT sdim,
|
||||
LADOU wr, LADOU wi, LADOU vs, CONST_LAINT ldvs,
|
||||
LADOU work, CONST_LAINT lwork, const void* bwork, LAINT info);
|
||||
void LAPACK_dgecon(LACHAR norm, CONST_LAINT n, CONST_LADOU a, CONST_LAINT lda,
|
||||
CONST_LADOU anorm, LADOU rnorm, LADOU work, LAINT iwork,
|
||||
LAINT info);
|
||||
void LAPACK_dtrexc(LACHAR compq, CONST_LAINT n, LADOU t, CONST_LAINT ldt,
|
||||
LADOU q, CONST_LAINT ldq, LAINT ifst, LAINT ilst, LADOU work,
|
||||
LAINT info);
|
||||
void LAPACK_dtrsyl(LACHAR trana, LACHAR tranb, CONST_LAINT isgn, CONST_LAINT m,
|
||||
CONST_LAINT n, CONST_LADOU a, CONST_LAINT lda, CONST_LADOU b,
|
||||
CONST_LAINT ldb, LADOU c, CONST_LAINT ldc, LADOU scale,
|
||||
LAINT info);
|
||||
void LAPACK_dpotrf(LACHAR uplo, CONST_LAINT n, LADOU a, CONST_LAINT lda,
|
||||
LAINT info);
|
||||
void LAPACK_dgges(LACHAR jobvsl, LACHAR jobvsr, LACHAR sort, DGGESCRIT delztg,
|
||||
CONST_LAINT n, LADOU a, CONST_LAINT lda, LADOU b, CONST_LAINT ldb,
|
||||
LAINT sdim, LADOU alphar, LADOU alphai, LADOU beta,
|
||||
LADOU vsl, CONST_LAINT ldvsl, LADOU vsr, CONST_LAINT ldvsr,
|
||||
LADOU work, CONST_LAINT lwork, LAINT bwork, LAINT info);
|
||||
void LAPACK_dsyev(LACHAR jobz, LACHAR uplo, CONST_LAINT n, LADOU a, CONST_LAINT lda,
|
||||
LADOU w, LADOU work, CONST_LAINT lwork, LAINT info);
|
||||
void LAPACK_dsyevr(LACHAR jobz, LACHAR range, LACHAR uplo, CONST_LAINT n, LADOU a,
|
||||
CONST_LAINT lda, LADOU lv, LADOU vu, CONST_LAINT il, CONST_LAINT iu,
|
||||
CONST_LADOU abstol, LAINT m, LADOU w, LADOU z, CONST_LAINT ldz,
|
||||
LAINT isuppz, LADOU work, CONST_LAINT lwork, LAINT iwork, CONST_LAINT liwork,
|
||||
LAINT info);
|
||||
};
|
||||
|
||||
|
||||
#endif /* CPPLAPACK_H */
|
||||
|
||||
|
||||
// Local Variables:
|
||||
// mode:C++
|
||||
// End:
|
||||
|
||||
Loading…
Reference in New Issue