[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