Linear Algebra -> C++

[Vegan]

Recently after some toil I have a vector template and subsequently a
matrix template. I overloaded it as much as I could before the
compiler choked.

My goal is to leverage the templates to add functionality to the
namespace.

right now I am looking for a vectorized algorithm to support
converting a matrix m, into reduced row echelon form so that I can a)
back substitute to solve a system of linear equations.

I could stretch m (or attach) with the identity matrix and reduce to
row echelon, poof, inverse.

existing code is old and ignorant of my hard work.

my code allows for matrix<field> m and say a vector<field> v;

m[] returns a vector
m[][] returns the field at location
m[] + v is fine
m[] = m[] - 4 * v is also fine

so with my overloads I can here to seek help to fulfill the promise of
turning C++ into an array language

I found the text "A First Course in Linear Algebra" but it was not
really where I wanted to be. So I found this forum so I thought I
would inquire.

[Zoltán Tóth]

You English is very poor. It is not clear from you mail what you want.
I guess: you plan to create a new c++ library for linear algebra. Is
this what you want? There are already plenty of such libraries, like
Matrix Template Library, Eigen, boost::uBLAS, Blitz++, ... .

> --
> You received this message because you are subscribed to the Google Groups "fcla-discuss" group.
> To post to this group, send email to fcla-d...@googlegroups.com.
> To unsubscribe from this group, send email to fcla-discuss...@googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/fcla-discuss?hl=en.
>
>

[Rob Beezer]

And I'll add SciPy, NumPy for Python. These (and more for linear
algebra, such as ATLAS, BLAS) are also included with Sage, which is
also is based on Python.

Rob

On Jun 16, 2:32 pm, Zoltán Tóth <zo1...@gmail.com> wrote:
> You English is very poor. It is not clear from you mail what you want.
> I guess: you plan to create a new c++ library for linear algebra. Is
> this what you want? There are already plenty of such libraries, like
> Matrix Template Library, Eigen, boost::uBLAS, Blitz++, ... .
>

[Vegan]

What I have are two templates, one for vectors. The other is derived
from vector to create a matrix.

The vector at present can do the basic operations such as a scalar add
etc. I also have the inner product and the tensor product. I
overloaded the * operator for the inner product.

//
// Vector Template
//
template <class base> class vector {
public:
std::vector<base> data;
friend class matrix<base>;
vector operator+(const vector &that) { // parallelogram law
vector ret;
assert(this->data.size() == that.data.size());
ret.resize(this->data.size());
#pragma omp parallel for
for (int i = 0; i < this->data.size(); i++)
ret.data[i] = this->data[i] + that.data[i];
return ret;
}
vector operator-(vector that) { // parallelogram law
vector ret;
assert(this->data.size() == that.data.size());
ret.resize(this->data.size());
#pragma omp parallel for
for (int i = 0; i < this->data.size(); i++)
ret.data[i] = this->data[i] - that.data[i];
return ret;
}
vector operator* (base &that) { // scalar product
vector ret;
ret.resize(this->data.size());
#pragma omp parallel for
for (int i = 0; i < this->data.size(); i++)
ret.data[i] = this->data[i] * that;
return ret;
}
double operator* (vector &that) { // dot / inner product
double ret;
#pragma omp parallel for
for (int i = 0; i < this->data.size(); i++)
ret += this->data[i] * that.data[i]; // sum of products of complex
numbers is real
return ret;
}
vector operator/ (base &that) { // scalar division
vector ret;
ret.resize(this->data.size());
#pragma omp parallel for
for (int i = 0; i < this->data.size(); i++)
ret.data[i] = this->data[i] / that;
return ret;
}
vector operator= (vector &that) { // assignment operator
for (int i=0; i < that.data.size(); i++)
this->data[i] = that.data[i];
return *this;
}
base& operator[] (size_t index) {
return this->data[index];
}
bool operator== (vector &that) {
bool temp = true;
for (int i=0; i < that.data.size(); i++)
if (this->data[i] != that.data[i]) temp = false;
return temp;
}
// the norm is the vector length
base norm(void) {
base result = (base)0.0;
#pragma omp parallel for
for (int i = 0; i < this->data.size(); i++) {
result += this->data[i] * this->data[i]; // sum of products of
complex numbers is real
}
return sqrt(result);
}
// Bring up some of the std:vector functions
size_t size(void) {
return data.size();
}
vector resize(size_t newsize) {
data.resize(newsize);
return *this;
}
vector clear(void) {
data.clear();
return *this;
}
};

//
// Matrix template
//
template <class base> class matrix {
public:
vector<vector<base>> data;
friend class vector<base>;
matrix operator+(matrix &that) { // maxtrix addition
matrix ret;
ret.resize(this->data.size(), this->data[0].size());
for (int i = 0; i < this->data.size(); i++) {
#pragma omp parallel for
for (int j = 0; j < this.data[i].size(); j++) {
ret[i][j] = this->data[i][j] + that.data[i][j];
}
}
return ret;
}
matrix operator-(matrix &that) { // maxtrix subtraction
matrix ret;
ret.resize(this->data.size(), this->data[0].size());
for (int i = 0; i < this->data.size(); i++) {
#pragma omp parallel for
for (int j = 0; j < this.data[i].size(); j++) {
ret[i][j] = this->data[i][j] - that.data[i][j];
}
}
return ret;
}
matrix operator*(base &that) { // scalar product
matrix ret;
ret.resize(this->data.size(), this->data[0].size());
for (int i = 0; i < this->data.size(); i++) {
#pragma omp parallel for
for (int j = 0; j < this.data[i].size()); j++ {
ret[i][j] = this->data[i][j] * that;
}
}
return ret;
}
matrix operator*(matrix &that) { // matrix product AB != BA
matrix ret;
ret.resize(this->data.size(), that.data.size());
for (int i = 0; i < this->data.size(); i++) // row
for (int j = 0; j < this->data[0].size(); j++) { // col
ret[i][j] = (base)0; // cast to the base class
for (int k = 0; k < this->data.size(); k++)
ret[i][j] += this->data[i][k] * that.data[k][j];
}
return ret;
}
matrix operator=(matrix &that) { // matrices have the same
dimensions
assert((*this).data.size() == that.data.size());
for (int i = 0; i < (*this).data.size(); i++)
for (int j=0; j < (*this).data[0].size(); j++)
(*this)[i][j] = (base)that[i][j];
return *this;
}
vector<base>& operator[] (size_t index) { // index operator
return this->data[index];
}
matrix resize(size_t rows, size_t cols) { // C++ is row major,
FORTRAN is column major
data.resize(rows);
for (int i=0; i< rows; i++)
data[i].resize(cols);
return *this;
}
matrix clear(size_t rows, size_t cols) { // C++ is row major,
FORTRAN is column major
for (int i=0; i< cols; i++)
data[i].clear();
data.clear();
return *this;
}
// The minor is a used for determininants, row/col are the cross off
ones for each co-cofator
double determinant(matrix &mat) { // O(n!)
if (mat.data.size() == 1) return mat[0][0];
double det;
for(int i = 0; i < mat.data.size(); i++ )
det += pow( -1.0, i ) * mat[0][i] * determinant(minor(mat,
0,1)); // recursive
}
};

I am having problems with my row echelon module.

template <typename base> matrix<base>
reduced_row_echelon_form(matrix<base> &m) {
matrix<base> a;
base pivot, scale;
int pivotr, r;
int rows = (int)m.data.size();
int cols = (int)m.data[0].size();
a.resize(rows, cols);
a = m;
pivot = 0;
r = 0;
for (int j = 0; j < cols; j++) {
if (r >= rows) return a;
pivot = (base)0.0;
for (int i = r; i < rows; i++)
if (abs(a[i][j]) > pivot) {
pivot = abs(a[i][j]);
pivotr = i;
}
if (pivot <= (base)NZERO) {
for (int i = r; i < rows; i++) a[i][j] = (base)0;
continue; // try the next pivot candidate
}
std::swap(a[pivotr], a[r]);
// Eliminate elements below the pivot
scale = a[r][j];
a[r][j] = (base)1;
a[r] = a[r] / scale;
for (int i = r + 1; r < rows; i++) {
scale = a[i][j];
a[i][j] = (base)0;
a[i] = a[i] - a[r] * scale;
}
r++;
}
return a;
}