99 views

Skip to first unread message

Jun 16, 2010, 1:29:08 PM6/16/10

to fcla-discuss

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.

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.

Jun 16, 2010, 5:32:49 PM6/16/10

to fcla-d...@googlegroups.com

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++, ... .

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.

>

>

Jun 16, 2010, 7:51:47 PM6/16/10

to fcla-discuss

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++, ... .

>

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++, ... .

>

Jun 16, 2010, 10:47:19 PM6/16/10

to fcla-discuss

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;

}

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;

}

Reply all

Reply to author

Forward

0 new messages

Search

Clear search

Close search

Google apps

Main menu