spence@psych.toronto.edu (Ian Spence) (11/21/89)
Has anyone had experience with either: (a) The CommonMatrix class library from ImageSoft, Inc (b) The M++ matrix class library from Ansys software. Any comparisons, or opinions regarding either, would be welcome. ---------------------------------------------------------------- Ian Spence spence@psych.utoronto.ca Department of Psychology spence@psych.toronto.edu University of Toronto SPENCE@UTORVM (BITNET) Toronto, Ontario (416) 978-7623 (Voice) Canada M5S 1A1 (416) 978-4811 (FAX)
sri@milton.u.washington.edu (Kandiah Sribalaskandarajah) (02/05/91)
Friends, are there any public domain matrix class libraries available? Thanks. ---
ericg@ucschu.ucsc.edu (Eric Goodman) (03/17/91)
Are there any publically available? If so, where? TIA (Thanks in advance) Eric Goodman, UC Santa Cruz ericg@ucschu.ucsc.edu or @ucschu.bitnet Eric_Goodman.staff@macmail.ucsc.edu
bjaspan@athena.mit.edu (Barr3y Jaspan) (03/18/91)
I've had several requests for this, and it isn't very long, so I guess I
should just post it here. I have found it very useful and reliable; your
results may vary. Naturally, no docs are provided (UTFL). :-) I don't claim
that is the most efficient or well designed or, in fact, anything.
Please send any extensions/bugfixes back to me.
--
Barr3y Jaspan, bjaspan@mit.edu
Watchmaker Computing
---- snip Matrix.H snip ----
#ifndef __MATRIX_H__
#define __MATRIX_H__
/* $Id: Matrix.H,v 1.1 90/12/03 03:23:20 marc Exp Locker: marc $ */
#include <assert.h> /* for abort() */
#include <stream.h>
class Matrix {
public:
/* Constructors */
Matrix(int, int);
Matrix(int, int, int *);
Matrix(int, int, double *);
Matrix(int, int, int ...);
Matrix(int, int, double ...);
Matrix(const Matrix&);
Matrix(const Matrix&,int,int,int,int);
~Matrix();
/* Member functions */
inline int Rows() const;
inline int Cols() const;
inline int SameSize(const Matrix&) const;
inline double& operator()(int, int) const;
inline double *Row(int) const;
void Print() const;
/* Overloaded operators */
friend Matrix operator+(const Matrix&, const Matrix&);
friend Matrix operator-(const Matrix&, const Matrix&);
friend Matrix operator*(const Matrix&, const Matrix&);
friend Matrix operator*(double, const Matrix&);
Matrix& operator+=(const Matrix&);
Matrix& operator-=(const Matrix&);
inline friend Matrix& operator*=(Matrix&, const Matrix&);
inline friend Matrix operator/(const Matrix&, double);
friend int operator==(const Matrix&, const Matrix&);
friend Matrix operator|(const Matrix&, const Matrix&);
friend Matrix operator&(const Matrix&, const Matrix&);
Matrix& operator=(const Matrix&);
friend ostream& operator<<(ostream&, const Matrix&);
/* Special matrix creators */
static inline const Matrix& NAM();
/* Matrix mutators */
void Identify();
private:
void Init(int, int);
static Matrix *M_NAM;
static Matrix *MakeNAM();
static double RowDotColumn(const Matrix&, int i, const Matrix&, int j);
int M_m, M_n, nam;
double *d;
};
inline int Matrix::Rows() const
{
return M_m;
}
inline int Matrix::Cols() const
{
return M_n;
}
inline Matrix& operator*=(Matrix& m1, const Matrix& m)
{
return(m1 = m1*m);
}
inline Matrix operator/(const Matrix& m, double k)
{
return((1.0/k)*m);
}
inline double& Matrix::operator()(int m, int n) const
{
/* My kingdom for exception handling ... */
if (m > M_m || n > M_n || m < 1 || n < 1)
abort();
return d[M_n*(m-1) + n - 1];
}
inline double *Matrix::Row(int r) const
{
return(d+(r-1)*M_n);
}
inline const Matrix& Matrix::NAM()
{
return *Matrix::M_NAM;
}
inline int Matrix::SameSize(const Matrix& m) const
{
return (M_m == m.M_m && M_n == m.M_n);
}
/* DO NOT ADD ANYTHING AFTER THIS #endif */
#endif /* __MATRIX_H__ */
---- snip Matrix.C snip ----
/*
* C++ class for MxN matrices.
*
* $Id: Matrix.C,v 1.2 90/12/04 01:38:15 bjaspan Exp Locker: marc $
*
*/
/* An m x n matrix has m rows and n columns -- ie height x width */
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include "Matrix.H"
Matrix *Matrix::M_NAM = Matrix::MakeNAM();
/* Constructors */
void Matrix::Init(int m, int n)
{
M_m = m;
M_n = n;
nam = 0;
d = new double[m*n];
}
Matrix::Matrix(int m, int n)
{
Init(m, n);
}
Matrix::Matrix(int m, int n, double first ...)
{
int i;
va_list ap;
Init(m, n);
d[0] = first;
va_start(ap, first);
for (i=1; i < m*n; i += 1)
d[i] = va_arg(ap, double);
}
Matrix::Matrix(int m, int n, int first ...)
{
int i;
va_list ap;
Init(m, n);
d[0] = double(first);
va_start(ap, first);
for (i=1; i < m*n; i += 1)
d[i] = double(va_arg(ap, int));
}
Matrix::Matrix(int m, int n, int *array)
{
int i;
Init(m, n);
for (i=0; i < (m*n); i += 1)
d[i] = double(array[i]);
}
Matrix::Matrix(int m, int n, double *array)
{
int i;
Init(m, n);
for (i=0; i < (m*n); i += 1)
d[i] = array[i];
}
Matrix::Matrix(const Matrix& m1)
{
int i;
Init(m1.M_m, m1.M_n);
for (i=0; i < (M_m*M_n); i += 1)
d[i] = m1.d[i];
}
Matrix::Matrix(const Matrix& m1, int r, int c, int m, int n)
{
int i,j;
if (r<1 || c<1 || m1.Rows() < r-1+m || m1.Cols() < c-1+n)
nam = 1;
else {
Init(m,n);
for (i=1;i<=m;i++)
for (j=1;j<=n;j++)
(*this)(i,j) = m1(r+i-1,c+j-1);
}
}
Matrix::~Matrix()
{
delete [M_m*M_n] d;
}
/* Overloaded operators */
int operator==(const Matrix& m1, const Matrix& m2)
{
int i;
if (m1.nam != m2.nam)
return 0;
if (m1.nam == 1)
return 1;
if (! m1.SameSize(m2))
return 0;
for (i = 0; i < m1.M_m*m1.M_n; i++)
if (m1.d[i] != m2.d[i])
return 0;
return 1;
}
Matrix operator+(const Matrix& m1, const Matrix& m2)
{
if (m1.M_m != m2.M_m || m1.M_n != m2.M_n ||
m1 == Matrix::NAM() || m2 == Matrix::NAM())
return Matrix::NAM();
Matrix m(m1);
for (int i=0; i < (m1.M_m*m1.M_n); i++)
m.d[i] += m2.d[i];
return m;
}
Matrix operator-(const Matrix& m1, const Matrix& m2)
{
if (m1.M_m != m2.M_m || m1.M_n != m2.M_n ||
m1 == Matrix::NAM() || m2 == Matrix::NAM())
return Matrix::NAM();
Matrix m(m1);
for (int i=0; i < (m1.M_m*m1.M_n); i++)
m.d[i] -= m2.d[i];
return m;
}
Matrix operator*(const Matrix& m1, const Matrix& m2)
{
if (m1.M_n != m2.M_m || m1 == Matrix::NAM() || m2 == Matrix::NAM())
return Matrix::NAM();
Matrix m(m1.M_m, m2.M_n);
for (int i=1; i <= m1.M_m; i += 1)
for (int j=1; j <= m2.M_n; j += 1)
for (int k=1; k <= m1.M_n; k += 1)
m(i, j) += m1(i, k) * m2(k, j);
return m;
}
Matrix& Matrix::operator+=(const Matrix& m)
{
if (M_m != m.M_m || M_n != m.M_n || *this == Matrix::NAM() ||
m == Matrix::NAM())
nam = 1;
else
for (int i=0; i < (M_m*M_n); i++)
d[i] += m.d[i];
return *this;
}
Matrix& Matrix::operator-=(const Matrix& m)
{
if (M_m != m.M_m || M_n != m.M_n || *this == Matrix::NAM() ||
m == Matrix::NAM())
nam = 1;
else
for (int i=0; i < (M_m*M_n); i++)
d[i] -= m.d[i];
return *this;
}
Matrix operator*(double k, const Matrix& m)
{
Matrix n(m);
int i;
for (i=0;i<n.M_m*n.M_n;i++) n.d[i]*=k;
return(n);
}
Matrix& Matrix::operator=(const Matrix& m)
{
if (M_m == m.M_m && M_n == m.M_n)
bcopy(m.d, d, M_m*M_n*sizeof(double));
else if (M_m*M_n == m.M_m*m.M_n) {
bcopy(m.d, d, M_m*M_n*sizeof(double));
M_m = m.M_m;
M_n = m.M_n;
} else {
delete [M_m*M_n] d;
M_m = m.M_m;
M_n = m.M_n;
d = new double[M_m*M_n];
bcopy(m.d, d, M_m*M_n*sizeof(double));
}
nam = m.nam;
return(*this);
}
Matrix operator|(const Matrix& m1, const Matrix& m2)
// Concatenates two matrices horizontally. Requires that m1 and m2
// have the same number of rows.
// [ A ] | [ B ] --> [ A | B ]
{
if (m1.Rows() != m2.Rows() || m1 == Matrix::NAM() || m2 == Matrix::NAM())
return Matrix::NAM();
Matrix m(m1.Rows(), m1.Cols() + m2.Cols());
int r, c;
for (r = 1; r <= m1.Rows(); r += 1)
for (c = 1; c <= m1.Cols(); c += 1)
m(r, c) = m1(r, c);
for (r = 1; r <= m1.Rows(); r += 1)
for (c = 1; c <= m2.Cols(); c += 1)
m(r, c + m1.Cols()) = m2(r, c);
return m;
}
Matrix operator&(const Matrix& m1, const Matrix& m2)
// Concatenates two matrices vertically. Requires that m1 and m2
// have the same number of columns.
// [ A ]
// [ A ] & [ B ] --> [ - ]
// [ B ]
{
if (m1.Cols() != m2.Cols() || m1 == Matrix::NAM() || m2 == Matrix::NAM())
return Matrix::NAM();
Matrix m(m1.Rows()+m2.Rows(), m1.Cols());
int r, c;
for (r = 1; r <= m1.Rows(); r += 1)
for (c = 1; c <= m1.Cols(); c += 1)
m(r, c) = m1(r, c);
for (r = 1; r <= m2.Rows(); r += 1)
for (c = 1; c <= m2.Cols(); c += 1)
m(r + m1.Rows(), c) = m2(r, c);
#if 0
// This doesn't work, and I don't know why.
bcopy(m1.d, m.d, m1.Rows()*m1.Cols()*sizeof(double));
bcopy(m2.d, m.d + m1.Rows()*m1.Cols()*sizeof(double),
m2.Rows()*m2.Cols()*sizeof(double));
#endif
return m;
}
ostream& operator<<(ostream& s, const Matrix& m)
{
int i,j;
if (m.nam) {
s << "NAM\n";
return(s);
}
for (i = 1; i <= m.M_m; i++) {
for (j = 1; j <= m.M_n; j++)
#ifdef __GNUG__
s.form("%6.2f ", m(i, j));
#else
s << form("%6.2f ", m(i, j));
#endif
s << "\n";
}
return(s);
}
/* Member functions */
void Matrix::Print() const
{
cout << *this;
}
/* Special matrix constructors */
Matrix *Matrix::MakeNAM()
{
Matrix *m = new Matrix(1,1);
m->nam = 1;
return m;
}
/* Matrix mutators */
void Matrix::Identify()
{
int i;
if (M_m != M_n) {
nam = 1;
return;
}
bzero((char *) d, sizeof(double)*M_m*M_n);
for (i = 1; i <= M_m; i++)
(*this)(i, i) = 1;
}Ari.Huttunen@hut.fi (Ari Juhani Huttunen) (03/19/91)
In article <1991Mar18.025539.18998@athena.mit.edu> bjaspan@athena.mit.edu (Barr3y Jaspan) writes: >class Matrix { ... > inline double& operator()(int, int) const; ... >}; I think this points out one problem in C++. You really should be able to say: inline double& operator[][](int, int) const; The problem isn't very serious, but it just isn't natural to access array members using notation array(x,y), but array[x][y]. -- Ari Huttunen, email: Ari.Huttunen@hut.fi, phone: 90-7285944
bjaspan@athena.mit.edu (Barr3y Jaspan) (03/19/91)
In article <ARI.HUTTUNEN.91Mar18173446@hydroman.hut.fi>, Ari.Huttunen@hut.fi (Ari Juhani Huttunen) writes: |> I think this points out one problem in C++. You really should be able to |> say: |> inline double& operator[][](int, int) const; |> |> The problem isn't very serious, but it just isn't natural to access array |> members using notation array(x,y), but array[x][y]. At the risk of starting a religous war, I'll respond. :-) In this particular case, the first method is preferable. Being able to say array[x][y] implies (in the conventional C style) that array[x] is meaningful (it always is in C). For my Matrix class, it isn't meaningful (you could argue that I should have derived Matrix from Vector, in which case array[x] could be meaningful, but that isn't the point.) So if you really always wanted to use [] for array indexing and () only for function calls, [] would have to accept any number and type of arguments (ie: operator[](int, int).) But then it would be semantically identical to () anyway, so what's the point? In any case, I agree that it isn't a serious problem. -- Barr3y Jaspan, bjaspan@mit.edu Watchmaker Computing
Ari.Huttunen@hut.fi (Ari Juhani Huttunen) (03/19/91)
In article <1991Mar18.213545.20966@athena.mit.edu> bjaspan@athena.mit.edu (Barr3y Jaspan) writes: >In article <ARI.HUTTUNEN.91Mar18173446@hydroman.hut.fi>, Ari.Huttunen@hut.fi (Ari Juhani Huttunen) writes: >|> I think this points out one problem in C++. You really should be able to >|> say: >|> inline double& operator[][](int, int) const; >|> >|> The problem isn't very serious, but it just isn't natural to access array >|> members using notation array(x,y), but array[x][y]. >In this particular case, the first method is preferable. Being able to say >array[x][y] implies (in the conventional C style) that array[x] is meaningful >(it always is in C). If array[] is not defined (array[][] is) then it certainly would be an error to call this undefined member function. So I don't see any problem here. Perhaps it is a religious issue after all. :( > For my Matrix class, it isn't meaningful (you could >argue that I should have derived Matrix from Vector, in which case array[x] >could be meaningful, but that isn't the point.) No, if you did that, there would be some performance penalties or at least I believe so. Correct me if I'm wrong. (If I am indeed wrong, this whole conversation is not needed.) >So if you really always wanted to use [] for array indexing and () only for >function calls, [] would have to accept any number and type of arguments (ie: >operator[](int, int).) But then it would be semantically identical to () >anyway, so what's the point? Well, why do we then have operator[] vs. operator() if if operator() alone would be sufficient? Because operator[] is more clear with arrays. -- Ari Huttunen, email: Ari.Huttunen@hut.fi, phone: 90-7285944
sking@nowhere.uucp (Steven King) (03/20/91)
In article <1991Mar18.213545.20966@athena.mit.edu> bjaspan@athena.mit.edu (Barr3y Jaspan) writes: >In article <ARI.HUTTUNEN.91Mar18173446@hydroman.hut.fi>, Ari.Huttunen@hut.fi (Ari Juhani Huttunen) writes: >|> I think this points out one problem in C++. You really should be able to >|> say: >|> inline double& operator[][](int, int) const; >|> >|> The problem isn't very serious, but it just isn't natural to access array >|> members using notation array(x,y), but array[x][y]. > >At the risk of starting a religous war, I'll respond. :-) > >In this particular case, the first method is preferable. Being able to say >array[x][y] implies (in the conventional C style) that array[x] is meaningful >(it always is in C). For my Matrix class, it isn't meaningful (you could >argue that I should have derived Matrix from Vector, in which case array[x] >could be meaningful, but that isn't the point.) > >So if you really always wanted to use [] for array indexing and () only for >function calls, [] would have to accept any number and type of arguments (ie: >operator[](int, int).) But then it would be semantically identical to () >anyway, so what's the point? > Okay, at the risk of fanning the flames of discussion... why not: class Matrix { double *d ; int M_m, M_n ; public: class aux { double *vec ; aux ( double *v ) { vec = v ; } friend class Matrix ; public: double& operator [] ( unsigned n ) { return vec [ n ] ; } } ; aux operator [] ( unsigned n ) { return aux ( &d [ n * M_m ] ) ; } } ; Obviously, various bounds checking etc., should be included. This then does allow for statements of the form: Matrix m ; m [j][i] = 1.0 ; Of course it does create an extra temporary and and a copy assignment of the matrix aux, but it does give that syntactical sugaring... >In any case, I agree that it isn't a serious problem. > >-- >Barr3y Jaspan, bjaspan@mit.edu >Watchmaker Computing -- Look Ma! No .sig! ..!cs.utexas.edu!ut-emx!nowhere!sking
stephens@motcid.UUCP (Kurt Stephens) (03/20/91)
Ari.Huttunen@hut.fi (Ari Juhani Huttunen) writes: >In article <1991Mar18.025539.18998@athena.mit.edu> bjaspan@athena.mit.edu (Barr3y Jaspan) writes: >>class Matrix { > ... >> inline double& operator()(int, int) const; > ... >>}; >I think this points out one problem in C++. You really should be able to >say: > inline double& operator[][](int, int) const; >The problem isn't very serious, but it just isn't natural to access array >members using notation array(x,y), but array[x][y]. Just define: inline double* Matrix::operator [](int i) { return Row(i); } Matrix M; double d; d = M[1][2]; "M[1]" evals to a double*, right? The "[2]" indexes off the double*. KISS (Keep It Simple Stephens) ;^) -- Kurt A. Stephens Foo::Foo(){return Foo();} stephens@void.rtsg.mot.com "When in doubt, recurse."
ahodgson@athena.mit.edu (Antony Hodgson) (03/20/91)
In article <ARI.HUTTUNEN.91Mar19102801@silver-surfer.hut.fi> Ari.Huttunen@hut.fi (Ari Juhani Huttunen) writes: >In article <1991Mar18.213545.20966@athena.mit.edu> bjaspan@athena.mit.edu (Barr3y Jaspan) writes: stuff about whether matrices should be indexed by (x,x) or [x][x]. > >If array[] is not defined (array[][] is) then it certainly would be an error >to call this undefined member function. So I don't see any problem here. >Perhaps it is a religious issue after all. :( > >> For my Matrix class, it isn't meaningful (you could >>argue that I should have derived Matrix from Vector, in which case array[x] >>could be meaningful, but that isn't the point.) > >No, if you did that, there would be some performance penalties or at least >I believe so. Correct me if I'm wrong. (If I am indeed wrong, this whole >conversation is not needed.) > My matrix class overloads [] to return a Vector&. Because of the way my code is structured, there's no performance penalty (bounds checking takes place, but that has to happen anyway). The only problem is that someone outside the class can change the size of the Vector returned. I'm not sure I can avoid that and still allow the elements of the vector to function as LHSides (i.e., Vector overloads [] to return a double& so that individual elements can be changed). I suppose one could derive a class ProtectedVector privately from vector and disallow size-changing operations on it; Matrix [] could then return a reference to a protected vector. Tony Hodgson ahodgson@hstbme.mit.edu
ahodgson@athena.mit.edu (Antony Hodgson) (03/21/91)
In article <6732@celery34.UUCP> stephens@motcid.UUCP (Kurt Stephens) writes: >Ari.Huttunen@hut.fi (Ari Juhani Huttunen) writes: stuff about how he wants to use [][] to get access to matrix elements rather than (int,int). > >Just define: > >inline >double* >Matrix::operator [](int i) { return Row(i); } > >Matrix M; >double d; > > d = M[1][2]; > >"M[1]" evals to a double*, right? The "[2]" indexes off the double*. >KISS (Keep It Simple Stephens) ;^) > In this case, the solution is too simple. The problem with returning a double* is that you can no longer perform bounds checking on the second index (the index i in operator[] should really be checked for validity as well). Of course, if you don't want to do bounds checking, this solution is fine. Tony Hodgson ahodgson@hstbme.mit.edu