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...
--
Opera: Sing it loud! :o( )>-<
/***************************************************************************
* Copyright (C) 2006 by Alan Grimes *
* [EMAIL PROTECTED] *
* *
* 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
/***************************************************************************
* Copyright (C) 2006 by Alan Grimes *
* [EMAIL PROTECTED] *
* *
* 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;
}
// ###################################
/***************************************************************************
* Copyright (C) 2006 by Alan Grimes *
* [EMAIL PROTECTED] *
* *
* 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
/***************************************************************************
* Copyright (C) 2006 by Alan Grimes *
* [EMAIL PROTECTED] *
* *
* 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@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/ktechlab-devel
