[Ktechlab-devel] Matrix and vector classes...

Jason Lucas jason.lucas7 at ntlworld.com
Wed Jun 27 03:11:02 UTC 2007


Alan,

I'll take a look, however I think it may be more helpful to the
community in general to add your fork to a branch in the repository.

Jason

On Mon, 2007-06-25 at 23:41 -0400, Alan Grimes wrote:
> Note, in Ktechlab there is a class Matrix *AND* a class matrix. =P
> 
> the issue I have with how they're done is that the matrix implementation
> is more than a bit messy and it uses a list of vectors internally,
> vectors are a fairly heavy-weight objects in themselves and this adds
> overhead.
> 
> This is a set of matrix and vector classes I develed for a class I was
> taking. I thought I could just drop them into ktechlab... Well, they
> were written using Smalltalk conventions while Ktechlab uses C++
> conventions... It took several days to do it....
> 
> Go ahead and make a "math" subdirectory for the project and put these in
> it. Remove the classes matrix and Vector and use these instead.
> 
> They have a number of methods which should be helpful in cleaning up the
> code of Matrix and elsewhere.
> 
> I have been using these in my own fork... This is part of the task of
> reconciling the changes...
> 
> plain text document attachment (qmatrix.h)
> /***************************************************************************
>  *   Copyright (C) 2006 by Alan Grimes   *
>  *   agrimes at speakeasy.net   *
>  *                                                                         *
>  *   This program is free software; you can redistribute it and/or modify  *
>  *   it under the terms of the GNU General Public License as published by  *
>  *   the Free Software Foundation; either version 2 of the License, or     *
>  *   (at your option) any later version.                                   *
>  *                                                                         *
>  *   This program is distributed in the hope that it will be useful,       *
>  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
>  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
>  *   GNU General Public License for more details.                          *
>  *                                                                         *
>  *   You should have received a copy of the GNU General Public License     *
>  *   along with this program; if not, write to the                         *
>  *   Free Software Foundation, Inc.,                                       *
>  *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
>  ***************************************************************************/
> 
> #ifndef QMATRIX_H
> #define QMATRIX_H
> 
> #ifndef QVECTOR_H
> #include "qvector.h"
> #endif
> 
> class QuickMatrix {
> public :
> 	 QuickMatrix(CUI m_in, CUI n_in);
> 	 QuickMatrix(CUI m_in);
> 	 QuickMatrix(const QuickMatrix *old); // ye olde copy constructor.
> 	~QuickMatrix();
> 
> // accessors
> // we use accessors so that we can provide range checking.
> // we use Smalltalk style naming conventions.
> 	double at(CUI m_a, CUI n_a) const;
> 	bool atPut(CUI m_a, CUI n_a, const double val);
> 	bool atAdd(CUI m_a, CUI n_a, const double val); // just give a negative val to subtract. =)
> 
> 	bool isSquare() const;
> 
> 	double *&operator[]( const int i) { return values[i]; }
> 	const double *operator[]( const int i) const { return values[i]; }
> 
> 	unsigned int size_m() const;
> 	unsigned int size_n() const;
> 
> // functions for some elementary row operations.
> // these are actually procedures because they operate on the current matrix rather than
> // producing a results matrix.  
> 	bool scaleRow(CUI m_a, const double scalor);
> 		// changes B by adding A.
> 	bool addRowToRow(CUI m_a, CUI m_b);
> 		// changes B by adding the result of A times a scalor 
> 	bool scaleAndAdd(CUI m_a, CUI m_b, const double scalor);
> 	bool partialScaleAndAdd(CUI m_a, CUI m_b, const double scalor);
> 	bool partialSAF(CUI m_a, CUI m_b, CUI from, const double scalor);
> 	bool swapRows(CUI m_a, CUI m_b);
> 
> // functions that accelerate certain types of
> // operations that would otherwise require millions of at()'s
> 	double multstep(CUI row, CUI pos, CUI col) const;
> 	double multRowCol(CUI row, CUI col, CUI lim) const;
> 
> 	QuickMatrix *transposeSquare() const; // Multiplies self by transpose.
> 	QuickVector *transposeMult(const QuickVector *operandvec) const;
> 
> // utility functions:
> 	void fillWithRandom();
> 	void fillWithZero();
> 	double rowsum(CUI m);
> 	double absrowsum(CUI m);
> 	QuickVector *normalizeRows();
> 
> // Matrix arithmetic.
> 	QuickMatrix *operator +=(const QuickMatrix *othermat);
> 	QuickMatrix *operator *=(const double y);
> 	QuickMatrix *operator =(const double y); // sets the diagonal to a constant.
> 	QuickMatrix *operator =(const QuickMatrix *oldmat);
> 	QuickVector *operator *(const QuickVector *operandvec) const;
> 	QuickMatrix *operator *(const QuickMatrix *operandmat) const;
> 
> // debugging
> 	void dumpToAux() const;
> 
> private :
> // We don't have a default matrix size so therefore we lock the default constructor. 
> 	QuickMatrix() {};
> 
> 	void allocator();
> 
> 	unsigned int m, n;
> 	double **values;
> };
> 
> #endif
> plain text document attachment (qvector.cpp)
> /***************************************************************************
>  *   Copyright (C) 2006 by Alan Grimes   *
>  *   agrimes at speakeasy.net   *
>  *                                                                         *
>  *   This program is free software; you can redistribute it and/or modify  *
>  *   it under the terms of the GNU General Public License as published by  *
>  *   the Free Software Foundation; either version 2 of the License, or     *
>  *   (at your option) any later version.                                   *
>  *                                                                         *
>  *   This program is distributed in the hope that it will be useful,       *
>  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
>  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
>  *   GNU General Public License for more details.                          *
>  *                                                                         *
>  *   You should have received a copy of the GNU General Public License     *
>  *   along with this program; if not, write to the                         *
>  *   Free Software Foundation, Inc.,                                       *
>  *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
>  ***************************************************************************/
> 
> #include <cstdlib> // for null 
> #include <cmath>
> #include <cassert>
> #include <kdebug.h>
> #include <iostream>
> 
> using namespace std;
> using namespace std;
> 
> #ifndef QVECTOR_H
> #include "qvector.h"
> #endif
> 
> /*
> #ifndef BADRNG_H
> #include "badrng.h"
> #endif
> */
> 
> // ######################################
> 
> unsigned int QuickVector::size() const {
> 	return m;
> }
> 
> // ######################################
> 
> double QuickVector::at(CUI m_a) const {
> //	if(!m_a || m_a > m) return NAN;
> 
> 	return values[m_a];
> }
> 
> // #####################################
> 
> bool QuickVector::atPut(CUI m_a, const double val) {
> 	if(m_a >= m) return false;
> 
> 	values[m_a] = val;
> 	changed = true;
> 	return true;
> }
> 
> // #####################################
> 
> bool QuickVector::atAdd(CUI m_a, const double val) {
> 	if(m_a >= m) return false;
> 
> 	values[m_a] += val;
> 	changed = true;
> 	return true;
> }
> 
> // #####################################
> 
> QuickVector::QuickVector(CUI m_in)
> 	: m(m_in), changed(true)
> {
> 	assert(m);
> 	values = new double[m];
> 	memset(values, 0, sizeof(double) *m);
> }
> 
> // #####################################
> 
> QuickVector::QuickVector(const QuickVector *old)
> 	: m(old->m), changed(old->changed) {
> 	assert(m);
> 	values = new double[m];
> 
> 	for(unsigned int j = 0; j < m; j++) 
> 		values[j] = old->values[j];
> }
> 
> // #####################################
> 
> QuickVector::~QuickVector() {
> 	delete values;
> }
> 
> // #####################################
> 
> /*
> void QuickVector::fillWithRandom() {
> 	for(unsigned int j = 0; j < m; j++)
> 		values[j] = drng();
> }
> */
> 
> // #####################################
> 
> void QuickVector::fillWithZeros() {
> 	memset(values, 0, m*sizeof(double));
> 	changed = true;
> }
> 
> // #####################################
> 
> QuickVector *QuickVector::operator-(const QuickVector *y) const {
> 	if(y->m != m) return NULL;
> 
> 	QuickVector *ret;
> 	ret = new QuickVector(m);
> 
> 	for(unsigned int i = 0; i < m; i++) ret->values[i] = values[i] - y->values[i];
> 
> 	return ret;
> }
> 
> // ####################################
> 
> void QuickVector::dumpToAux() const {
> 	for(unsigned int i = 0; i < m; i++) cout << values[i] << ' ';
> 	cout << endl;
> }
> 
> // ####################################
> 
> bool QuickVector::swapElements(CUI m_a, CUI m_b) {
> 	if(m_a >= m || m_b >= m) return false;
> 
> 	double temp = values[m_a];
> 	values[m_a] = values[m_b];
> 	values[m_b] = temp;
> 	changed = true;
> 
> 	return true;
> }
> 
> // ###################################
> 
> QuickVector *QuickVector::operator *=(const QuickVector *y) {
> 	if(y->m != m) return NULL;
> 
> 	for(unsigned int i = 0; i < m; i++) values[i] *= y->values[i];
> 	changed = true;
> 	return this;
> }
> 
> // ###################################
> 
> QuickVector *QuickVector::operator +=(const QuickVector *y) {
> 	if(y->m != m) return NULL;
> 
> 	for(unsigned int i = 0; i < m; i++) values[i] += y->values[i];
> 	changed = true;
> 	return this;
> }
> 
> // ###################################
> 
> plain text document attachment (qvector.h)
> /***************************************************************************
>  *   Copyright (C) 2006 by Alan Grimes   *
>  *   agrimes at speakeasy.net   *
>  *                                                                         *
>  *   This program is free software; you can redistribute it and/or modify  *
>  *   it under the terms of the GNU General Public License as published by  *
>  *   the Free Software Foundation; either version 2 of the License, or     *
>  *(at your option) any later version.                                   *
>  *                                                                         *
>  *   This program is distributed in the hope that it will be useful,       *
>  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
>  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
>  *   GNU General Public License for more details.                          *
>  *                                                                         *
>  *   You should have received a copy of the GNU General Public License     *
>  *   along with this program; if not, write to the                         *
>  *   Free Software Foundation, Inc.,                                       *
>  *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
>  ***************************************************************************/
> 
> #ifndef QVECTOR_H
> #define QVECTOR_H
> 
> #ifndef CUI
> #define CUI	const unsigned int
> #endif
> 
> #define EPSILON	0.000001
> 
> class QuickVector {
> public :
> 	 QuickVector(CUI m_in);
> 	~QuickVector();
> 	 QuickVector(const QuickVector *old); // ye olde copy constructor.
> 
> 	double & operator[]( const int i) { changed = true; return values[i]; }
> 	const double operator[]( const int i) const { return values[i]; }
> 
> // accessors
> // we use accessors so that we can provide range checking.
> // we use Smalltalk style naming conventions.
> 	double at(CUI m_a) const;
> 	bool atPut(CUI m_a, const double val);
> 	bool atAdd(CUI m_a, const double val);
> 
> 	unsigned int size() const;
> 
> // utility functions:
> //	void fillWithRandom();
> 	void fillWithZeros();
> 	bool swapElements(CUI m_a, CUI m_b);
> 
> // Vector arithmetic.
> 
> 	QuickVector *operator*=(const double *y);
> 	QuickVector *operator*=(const QuickVector *y);
> 	QuickVector *operator+=(const QuickVector *y);
> 	QuickVector *operator-(const QuickVector *y) const;
> 
> // debugging
> 	void dumpToAux() const;
> 
> 	/**
> 	 * Returns true if the vector has changed since setUnchanged was last called
> 	 */
> 	inline bool isChanged() const { return changed; }
> 	/**
> 	 * Sets the changed status to false.
> 	 */
> 	inline void setUnchanged() { changed=false; }
> 
> private :
> // We don't have a default vector size so therefore we lock the default constructor.
> 	QuickVector() {};
> 	unsigned int m;
> 	bool changed;
> 	double *values;
> };
> 
> #endif
> plain text document attachment (qmatrix.cpp)
> /***************************************************************************
>  *   Copyright (C) 2006 by Alan Grimes   *
>  *   agrimes at speakeasy.net   *
>  *                                                                         *
>  *   This program is free software; you can redistribute it and/or modify  *
>  *   it under the terms of the GNU General Public License as published by  *
>  *   the Free Software Foundation; either version 2 of the License, or     *
>  *   (at your option) any later version.                                   *
>  *                                                                         *
>  *   This program is distributed in the hope that it will be useful,       *
>  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
>  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
>  *   GNU General Public License for more details.                          *
>  *                                                                         *
>  *   You should have received a copy of the GNU General Public License     *
>  *   along with this program; if not, write to the                         *
>  *   Free Software Foundation, Inc.,                                       *
>  *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
>  ***************************************************************************/
> 
> #include <cstdlib> // for NULL
> #include <cmath>
> #include <cassert>
> #include <iostream>
> 
> using namespace std;
> 
> #ifndef QMATRIX_H
> #include "qmatrix.h"
> #endif
> 
> /*
> #ifndef BADRNG_H
> #include "badrng.h"
> #endif
> */
> 
> // ####################################
> 
> bool QuickMatrix::isSquare() const {
> 	return m == n;
> } 
> 
> // ####################################
> 
> // Ideally, this should be inlined, however I am reluctant to put code in the header files which may be included in numerous
> // object files locations.
> 
> unsigned int QuickMatrix::size_m() const {
> 	return m;
> }
> 
> // ####################################
> 
> unsigned int QuickMatrix::size_n() const {
> 	return n;
> }
> 
> // ####################################
> 
> void QuickMatrix::allocator() {
> //	assert(!values);
> 	assert(m);
> 	assert(n);
> 
> 	values = new double*[m];
> 	for(unsigned int i = 0; i < m; i++) {
> 		values[i] = new double[n];
> 	}
> }
> 
> // ####################################
> 
> QuickMatrix *QuickMatrix::transposeSquare() const {
> 	QuickMatrix *newmat = new QuickMatrix(n);
> 
> 	for(unsigned int i= 0; i < n; i++) {
> 		for(unsigned int j = 0; j < n; j++) {
> 			newmat->values[i][j] = 0;
> 			for(unsigned int k = 0; k < m; k++) {
> 				newmat->values[i][j] += values[k][i] * values[k][j];
> 			}
> 		}
> 	}
> 
> 	return newmat;
> }
> 
> // #####################################
> 
> QuickVector *QuickMatrix::transposeMult(const QuickVector *operandvec) const {
> 	if(operandvec->size() != m) return NULL;
> 
> 	QuickVector *ret = new QuickVector(n);
> 
> 	for(unsigned int i = 0; i < n; i++) {
> 		double sum = 0;
> 		for(unsigned int j = 0; j < m; j++)
> 			sum += values[j][i] * (*operandvec)[j];
> 		(*ret)[i] = sum;
> 	}
> 
> 	return ret;
> }
> 
> // ####################################
> 
> QuickMatrix::QuickMatrix(CUI m_in, CUI n_in)
> 	: m(m_in), n(n_in) {
> 	allocator();
> 	fillWithZero();
> }
> 
> // ####################################
> 
> QuickMatrix::QuickMatrix(CUI m_in)
> 	: m(m_in), n(m_in) {
> 	allocator();
> 	fillWithZero();
> }
> 
> // ####################################
> 
> QuickMatrix::QuickMatrix(const QuickMatrix *old)
> 	: m(old->m), n(old->n) {
> 	allocator();
> 
> 	for(unsigned int j = 0; j < m; j++) {
> 		memcpy(values[j], old->values[j], n*sizeof(double)); // fastest method. =)
> 	}
> }
> 
> // ####################################
> 
> QuickMatrix::~QuickMatrix() {
> 	for(unsigned int i = 0; i < m; i++) delete values[i];
> 	delete values;
> }
> 
> // ####################################
> 
> double QuickMatrix::at(CUI m_a, CUI n_a) const
> {
> 	if(m_a >= m || n_a >= n) return NAN;
> 
> 	return values[m_a][n_a];
> }
> 
> // ####################################
> 
> double QuickMatrix::multstep(CUI row, CUI pos, CUI col) const
> {
> 	return values[row][pos] * values[pos][col];
> }
> 
> // ####################################
> 
> double QuickMatrix::multRowCol(CUI row, CUI col, CUI lim) const
> {
> 	const double *rowVals = values[row];
> 
> 	double sum = 0;
> 	for(unsigned int i = 0; i < lim; i++)
> 		sum += rowVals[i] * values[i][col];
> 	return sum;
> }
> 
> // ####################################
> 
> bool QuickMatrix::atPut(CUI m_a, CUI n_a, const double val) {
> 	if(m_a >= m || n_a >= n) return false;
> 
> 	values[m_a][n_a] = val;
> 	return true;
> }
> 
> // ####################################
> 
> bool QuickMatrix::atAdd(CUI m_a, CUI n_a, const double val) {
> 	if(m_a >= m || n_a >= n) return false;
> 
> 	values[m_a][n_a] += val;
> 	return true;
> }
> 
> // ####################################
> 
> bool QuickMatrix::swapRows(CUI m_a, CUI m_b) {
> 	if(m_a >= m || m_b >= m) return false;
> 
> 	double *temp = values[m_a];
> 	values[m_a] = values[m_b];
> 	values[m_b] = temp;
> 
> 	return true;
> }
> 
> // ####################################
> 
> bool QuickMatrix::scaleRow(CUI m_a, const double scalor) {
> 	if(m_a >= m) return false;
> 
> 	double *arow = values[m_a];
> 
> // iterate over n columns. 
> 	for(unsigned int j = 0; j < n; j++) arow[j] *= scalor;
> 
> 	return true;
> }
> 
> // ####################################
> 
> double QuickMatrix::rowsum(CUI m) {
> 	if(m >= m) return NAN;
> 
> 	double *arow = values[m];
> 
> 	double sum = 0.0;
> 
> // iterate over n columns. 
> 	for(unsigned int j = 0; j < n; j++) sum += arow[j];
> 
> 	return sum;
> }
> 
> // ####################################
> 
> double QuickMatrix::absrowsum(CUI m) {
> 	if(m >= m) return NAN;
> 
> 	double *arow = values[m];
> 
> 	double sum = 0.0;
> 
> // iterate over n columns. 
> 	for(unsigned int j = 0; j < n; j++) sum += fabs(arow[j]);
> 
> 	return sum;
> }
> 
> // ####################################
> 
> // behaves oddly but doesn't crash if m_a == m_b
> bool QuickMatrix::scaleAndAdd(CUI m_a, CUI m_b, const double scalor) {
> 	if(m_a >= m || m_b >= m) return false;
> 
> 	const double *arow = values[m_a];
> 	double *brow = values[m_b];
> 
> // iterate over n columns. 
> 	for(unsigned int j = 0; j < n; j++)
> 		brow[j] += arow[j] * scalor;
> 
> 	return true;
> }
> 
> // ####################################
> 
> // behaves oddly but doesn't crash if m_a == m_b
> bool QuickMatrix::partialScaleAndAdd(CUI m_a, CUI m_b, const double scalor) {
> 	if(m_a >= m || m_b >= m) return false;
> 
> 	const double *arow = values[m_a];
> 	double *brow = values[m_b];
> 
> // iterate over n - m_a columns.
> 	for(unsigned int j = m_a; j < n; j++)
> 		brow[j] += arow[j] * scalor;
> 
> 	return true;
> }
> 
> // ########################################
> 
> bool QuickMatrix::partialSAF(CUI m_a, CUI m_b, CUI from, const double scalor) {
> 	if(m_a >= m || m_b >= m) return false;
> 
> 	const double *arow = values[m_a];
> 	double *brow = values[m_b];
> 
> // iterate over n - m_a columns.
> 	for(unsigned int j = from; j < n; j++)
> 		brow[j] += arow[j] * scalor;
> 
> 	return true;
> }
> 
> // ####################################
> /*
> void QuickMatrix::fillWithRandom() {
> 	for(unsigned int j = 0; j < m; j++) {
> 		double *row = values[j];
> 		for(unsigned int i = 0; i < n; i++) row[i] = drng();
> 	}
> }
> */
> 
> // ####################################
> 
> QuickVector *QuickMatrix::normalizeRows() {
> 	QuickVector *ret = new QuickVector(m);
> 
> 	for(unsigned int j = 0; j < m; j++) {
> 		double *row = values[j];
> 		unsigned int max_loc = 0;
> 
> 		for(unsigned int i = 0; i < n; i++) {
> 			if(fabs(row[max_loc]) < fabs(row[i])) max_loc = i;
> 		}
> 
> 		double scalar = 1.0 / row[max_loc];
> 
> 		(*ret)[j] = scalar;
> 
> 		for(unsigned int i = 0; i < n; i++) row[i] *= scalar;
> 	}
> 
> 	return ret;
> }
> 
> // ####################################
> 
> QuickVector *QuickMatrix::operator *(const QuickVector *operandvec) const {
> 	if(operandvec->size() != n) return NULL;
> 
> 	QuickVector *ret = new QuickVector(m);
> 
> 	for(unsigned int i = 0; i < m; i++) {
> 		double sum = 0;
> 		for(unsigned int j = 0; j < n; j++)
> 			sum += values[i][j] * (*operandvec)[j];
> 		(*ret)[i] = sum;
> 	}
> 
> 	return ret;
> }
> 
> // ####################################
> 
> QuickMatrix *QuickMatrix::operator +=(const QuickMatrix *operandmat) {
> 	if(operandmat->n != n || operandmat->m != m) return NULL;
> 
> 	for(unsigned int i = 0; i < m; i++) {
> 		for(unsigned int j = 0; j < n; j++)
> 			values[i][j] += operandmat->values[i][j];
> 	}
> 
> 	return this;
> }
> 
> // ####################################
> 
> QuickMatrix *QuickMatrix::operator *(const QuickMatrix *operandmat) const {
> 	if(operandmat->m != n) return NULL;
> 
> 	QuickMatrix *ret = new QuickMatrix(m,operandmat->n);
> 
> 	for(unsigned int i = 0; i < m; i++) {
> 		for(unsigned int j = 0; j < operandmat->n; j++) {
> 			double sum = 0;
> 			for(unsigned int k = 0; k < n; k++) 
> 				sum += values[i][k] * operandmat->values[k][j];
> 			ret->values[i][j] = sum;
> 		}
> 	}
> 
> 	return ret;
> }
> 
> // ###################################
> 
> void QuickMatrix::dumpToAux() const {
> 	for(unsigned int j = 0; j < m; j++) {
> 		for(unsigned int i = 0; i < n; i++)
> 			cout << values[j][i] << ' ';
> 		cout << endl;
> 	} 
> }
> 
> // ###################################
> 
> void QuickMatrix::fillWithZero() {
> 	for(unsigned int j = 0; j < m; j++) {
> 		memset(values[j], 0, n*sizeof(double)); // fastest method. =)
> 	}
> }
> 
> // ###################################
> // sets the diagonal to a constant.
> QuickMatrix *QuickMatrix::operator =(const double y) {
> 	fillWithZero();
> 	unsigned int size = n;
> 	if(size > m) size = m;
> 
> 	for(unsigned int i = 0; i < size; i++) values[i][i] = y;
> 	return this;
> }
> 
> -------------------------------------------------------------------------
> This SF.net email is sponsored by DB2 Express
> Download DB2 Express C - the FREE version of DB2 express and take
> control of your XML. No limits. Just data. Click to get it now.
> http://sourceforge.net/powerbar/db2/
> _______________________________________________ Ktechlab-devel mailing list Ktechlab-devel at lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/ktechlab-devel





More information about the Ktechlab-devel mailing list