algebra

GenericMatrix.cpp

/*
 * GenericMatrix.cpp
 *
 *  Created on: Nov 20, 2015
 *      Author: lurker
 */

#include "GenericMatrix.h"

namespace outverse {

GenericMatrix::~GenericMatrix() {}

DenseMatrix::~DenseMatrix() {
	if(!_shared) mkl_free(_data);
}

void
DenseMatrix::zero() noexcept{
	memset(_data, 0., sizeof(double) * _m * _n);
}

void
DenseMatrix::multiply (GenericVector& v, GenericVector &r) const noexcept {
	cblas_dgemv(CblasRowMajor, CblasNoTrans, _m, _n, 1.0, _data, _m , v.get(), 1, 0.0, r.get(), 1);
}

void
DenseMatrix::multiply_t (GenericVector& v, GenericVector &r) const noexcept{
	cblas_dgemv(CblasRowMajor, CblasTrans,   _m, _n , 1.0, _data, _m, v.get(), 1, 0.0, r.get(), 1);
}

void
DenseMatrix::multiply(DenseMatrix &ma, DenseMatrix &rt) const noexcept {

}

void
DenseMatrix::multiply_t(DenseMatrix &ma, DenseMatrix &r) const noexcept {

}

GenericVector DenseMatrix::operator *(GenericVector& v) noexcept {
	GenericVector r(_m, true);
	if (!_transposed) {
		this->multiply(v, r);
	}
	else {
		this->multiply_t(v, r);
	}

	return r;
}

} /* namespace outverse */

GenericMatrix.h

/*
 * GenericMatrix.h
 *
 *  Created on: Nov 20, 2015
 *      Author: lurker
 */

#ifndef GENERICMATRIX_H_
#define GENERICMATRIX_H_


#include <cstring>
#include <cstdlib>
#include <vector>
#include <algorithm>
#include <cmath>
#include <mkl.h>
#include <cstdio>
#include <utility>
#include <cassert>
#include "GenericVector.h"

using namespace std;

namespace outverse {

class GenericMatrix {

friend class GenericVector;

public:
	GenericMatrix(size_t m, size_t n) noexcept {
		_m = m;
		_n = n;
	}
	virtual ~GenericMatrix() noexcept;

	size_t rows() noexcept {return _m;}
	size_t cols() noexcept {return _n;}

	virtual void multiply (GenericVector& v, GenericVector &r) const noexcept = 0;
	virtual void multiply_t (GenericVector& v, GenericVector &r) const noexcept = 0;
	virtual void zero() noexcept = 0;

protected:
	size_t _m;
	size_t _n;
};

class DenseMatrix : public GenericMatrix {

friend class GenericVector;

public:
	DenseMatrix(size_t m, size_t n) noexcept : GenericMatrix(m , n) {
		_shared = false;
		_transposed = false;
		_data =  (double *)mkl_malloc( _m * _n * sizeof( double ), 64 );
	}

	DenseMatrix(size_t m, size_t n, double* content) noexcept : GenericMatrix(m , n) {
		_shared = true;
		_transposed = false;
		_data = content;
	}

	virtual ~DenseMatrix();

	void multiply (GenericVector& v, GenericVector &r) const noexcept;
	void multiply_t (GenericVector& v, GenericVector &r) const noexcept;
	void multiply (DenseMatrix& ma, DenseMatrix & rt) const noexcept;
	void multiply_t (DenseMatrix& ma, DenseMatrix &rt) const noexcept;

	GenericVector operator* (GenericVector& v) noexcept;
	DenseMatrix operator* (DenseMatrix& ma) noexcept;

	double& operator() (size_t r, size_t c) noexcept {
		return *(_data + r * _n + c);
	}

	const double& operator() (size_t r, size_t c) const noexcept {
		return *(_data + r * _n + c);
	}

	void zero() noexcept;

	DenseMatrix& T() noexcept {
		_transposed = (!_transposed);
		return *this;
	}

private:
	double *_data;
	bool _shared;
	bool _transposed;
};

} /* namespace outverse */

#endif /* GENERICMATRIX_H_ */

GenericVector.cpp

/*
 * GenericVector.cpp
 *
 *  Created on: Nov 20, 2015
 *      Author: lurker
 */

#include "GenericVector.h"

namespace outverse {

GenericVector::GenericVector(size_t size, bool alloc) noexcept {
	_shared = false;
	_size   = size;

	if (alloc) {
		_data = (double *)mkl_malloc( _size * sizeof( double ), 64 );
	}
	else {
		_shared = true;
		_data = nullptr;
	}

}

GenericVector::GenericVector(size_t size, double *content, bool shared) noexcept {
	_size = size;
	if (shared) {
		_data = content;
		_shared = true;
	}
	else {
		_shared = false;
		_data   = (double *)mkl_malloc( _size * sizeof( double ), 64 );
		memcpy(_data, content, sizeof(double) * _size);
	}
}

const double &
GenericVector::operator [](size_t index) const noexcept {
	return *(_data + index);
}

double &
GenericVector::operator [](size_t index) noexcept {
	return *(_data + index);
}

bool
GenericVector::operator ==(const GenericVector& v) noexcept {
	if (v.size() != _size) return false;
	for (size_t i = 0; i < _size; ++i) {
		if (_data[i] != v[i]) {
			return false;
		}
	}
	return true;
}

size_t
GenericVector::size() const noexcept {
	return _size;
}

double
GenericVector::dot(const GenericVector& v) const noexcept {
	assert(v.size() == _size);
	return cblas_ddot(_size, _data, 1, v._data, 1);
}


double
GenericVector::dot(size_t size, const double* v) const noexcept {
	assert(size == _size);
	return cblas_ddot(_size, _data, 1, v , 1);
}

void
GenericVector::zero() noexcept {
	memset(_data, 0., sizeof(double) * _size);
}

double
GenericVector::norm() const noexcept {
	return cblas_dnrm2(_size, _data, 1);
}

void
GenericVector::normalize() noexcept {
	this->div(this->norm());
}

void
GenericVector::add(const GenericVector& v) noexcept {
	cblas_daxpy(_size, 1.0, v._data, 1, _data, 1);
}

void
GenericVector::sub(const GenericVector& v) noexcept {
	cblas_daxpy(_size, -1.0, v._data, 1, _data, 1);
}

void
GenericVector::mul(const double scalar) noexcept {
	cblas_dscal(_size, scalar, _data, 1);
}

void
GenericVector::div(const double scalar) noexcept {
	cblas_dscal(_size, 1.0/scalar, _data, 1);
}

void
GenericVector::operator +=(const GenericVector& v) noexcept {
	this->add(v);
}

void
GenericVector::operator -=(const GenericVector& v) noexcept {
	this->sub(v);
}

void
GenericVector::operator *=(double scalar) noexcept {
	this->mul(scalar);
}

void
GenericVector::operator /=(double scalar) noexcept {
	this->div(scalar);
}

void
GenericVector::assign(const GenericVector& v) noexcept {
	memcpy(_data, v._data, _size * sizeof(double));
}

GenericVector::~GenericVector() {
	if (!_shared) {
		mkl_free(_data);
	}
}

} /* namespace outverse */

GenericVector.h

/*
 * GenericVector.h
 *
 *  Created on: Nov 20, 2015
 *      Author: lurker
 */

#ifndef GENERICVECTOR_H_
#define GENERICVECTOR_H_

#include <cstring>
#include <cstdlib>
#include <vector>
#include <algorithm>
#include <cmath>
#include <mkl.h>
#include <cstdio>

#include <cassert>

using namespace std;

namespace outverse {

class GenericVector {
public:
	GenericVector(size_t size, bool alloc=true) noexcept;
	GenericVector(size_t size, double *content, bool shared=true) noexcept;
	const double& operator[](size_t index) const noexcept;
	double& operator[] (size_t index) noexcept;
	bool operator==(const GenericVector& v) noexcept;
	size_t size() const noexcept;

	double dot(size_t size, const double * const v) const noexcept;
	double dot(const GenericVector& v) const noexcept;
	void zero() noexcept;
	double norm() const noexcept;
	void normalize() noexcept;

	void add(const GenericVector& v) noexcept;
	void sub(const GenericVector& v) noexcept;
	void mul(const double scalar) noexcept;
	void div(const double scalar) noexcept;

	void operator+= (const GenericVector& v) noexcept;
	void operator-= (const GenericVector& v) noexcept;
	void operator*= (const double scalar) noexcept;
	void operator/= (const double scalar) noexcept;

	void assign(const GenericVector& v) noexcept;
	double* get() {return _data;}

	virtual ~GenericVector() noexcept;

private:
	double *_data;
	bool   _shared;
	size_t _size;
};

} /* namespace outverse */

#endif /* GENERICVECTOR_H_ */