|
Re: Need help - complex matrix, complex vector
Here it is. If you use the determinant function, you'll probably want to modify it to row-reduce and then multiply the diagonals. It uses the method of minors right now, because I was using matrices of polynomials, and so I couldn't do the division that is necessary to use the row-reduction method.
You can use this for vectors too. Just create a matrix with only one column. If you implement vectors that way, you can do matrix-vector multiplication.
If you distribute the code at all I would appreciate credit in the source. A link to my website ( www.blake-foster.com) would be nice.
EDIT: One more important note--all the indexes start from 1 instead of the usual 0. I did it that way to make it consistent with the notation used in math literature.
matrix.h:
#ifndef MATRIX_H
#define MATRIX_H
#include <vector>
#include <iostream>
#include <math>
using namespace std;
// The class Matrix implements matrices with entries in an arbitrary field.
// Rows, columns, or individual elements may be accessed in constant time. Rows
// and columns may be inserted only onto the end, and columns may be deleted
// from the end. Elements may be indexed in row major order as well.
// ***All indexes for matrices start with 1.***
template <typename T>
class Matrix
{
public:
Matrix<T>(int,int, T); // create int x int matrix filled T
Matrix<T>(int,int); // create int x int matrix filled with 0's
Matrix<T>(); // create an empty matrix
int numRows() const; // return the number of rows
int numCols() const; // return the number of columns
inline T* elem(int, int); // return a pointer to element (int,int).
Matrix<T> getRow(int); // return row int as a matrix
Matrix<T> getCol(int); // return column int as a matrix
void setRow(int, Matrix<T>); // takes a matrix with numCols elements and stores
//it in the (int)th column of THIS matrix.
void setCol(int, Matrix<T>); // takes a matrix with numRows elements and
//stores it in the (int)th row of THIS matrix.
Matrix<T> operator*(const Matrix<T>&); // matrix multiplication.
Matrix<T> operator*(const T&) const; // scalar multiplication.
Matrix<T> operator/(const T&) const; // inverse scalar multiplication.
// *** requires / operator to be implemented on T ***
bool operator==(const Matrix<T>&); // element by element equality
T* operator[](const int&); // returns the element of index int
// (upper left to lower right, row major order)
int length(); // the total number of elements
Matrix<T> transpose(); // the transpose of THIS matrix (swap rows and columns)
T determinant(); // the determinant
void push_row_back(T); // push a row filled with zeros
void push_col_back(T); // push a column filled with zeros
void pop_col_back(); // delete the last collumn
void pop_row_back(); // delete the last row
T cofactor(const int&, const int&); // compute the (i,j)th cofactor
Matrix<T> cofactor_matrix(); // compute cofactor matrix
Matrix<T> adjoint_matrix(); // compute adjoint matrix
Matrix<T> inverse(); // inverse matrix
// *** requires / operator to be implemented on T ***
void print();
protected:
int rows; // the number of rows
int cols; // the number of columns
vector<vector<T> > data; // the contents of the Matrix
};
template <typename T>
Matrix<T>::Matrix(int r, int c, T t) // construct an r x c matrix containing all zero's
: rows(r), cols(c), data(vector<vector<T> >(0))
{
for (int i=0; i < rows; ++i) data.push_back(vector<T>(cols,t));
}
template <typename T>
Matrix<T>::Matrix(int r, int c) // construct an r x c matrix containing all zero's
: rows(r), cols(c), data(vector<vector<T> >(0))
{
// initialize each row with a vector containing all zero's
for (int i=0; i < rows; ++i) data.push_back(vector<T>(cols,T()));
}
template <typename T>
Matrix<T>::Matrix() // construct an empty Matrix
{
Matrix<T>(0,0);
}
template <typename T>
int Matrix<T>::numRows() const {return rows;} // return the number of rows
template <typename T>
int Matrix<T>::numCols() const {return cols;} // return the number of columns
template <typename T>
inline T* Matrix<T>::elem(int r, int c) // return a pointer to the element in row r and
// column c (indexes start from 1)
{
return &data[r-1][c-1];
}
template <typename T>
Matrix<T> Matrix<T>::getRow(int index) // return row (index) as a matrix
{
Matrix<T> temp=Matrix<T>(1,cols); // creat a 1 x cols matrix
temp.data[0]=data[index-1]; // copy row (index) of THIS matrix into temp
return temp;
}
template <typename T>
Matrix<T> Matrix<T>::getCol(int index) // return column (index) as a matrix
{
Matrix<T> temp=Matrix<T>(rows,1); // creat a cols x 1 matrix
for (int i=1; i<=rows; ++i)
{
// copy elements of column (index) into temp one at a time
*temp.elem(i,1)=*elem(i,index);
}
return temp;
}
template <typename T>
void Matrix<T>::setRow(int index, Matrix<T> m) // replace row (index) of THIS Matrix
{
if (m.length()!=cols)
{
cout<<endl<<"dimension mismatch"<<endl; // m must have the same number of
// columns as THIS
system("pause");
return;
}
data[index-1]=m.data[0]; // copy m into row (index) of THIS Matrix
}
template <typename T>
void Matrix<T>::setCol(int index, Matrix<T> m)
{
if (m.length()!=rows)
{
cout<<endl<<"dimension mismatch"<<endl; // m must have the same number of
// rows as THIS
system("pause");
return;
}
for (int i=1; i<=rows; ++i) // iterate over all rows
{
*elem(i,index)=*m[i]; // copy the ith element of m into the ith position in
// row (index) of THIS
}
}
template <typename T>
Matrix<T> Matrix<T>::operator*(const Matrix<T>& m) // Matrix multiplication
{
Matrix b=m; // make a copy of m
if (cols != b.rows) // To multiply two matrices the number of colums of the
// first must match the number of columns of the second
{
cout << "dimension mismatch" << endl;
system("pause");
return *this;
}
T acc=T()-T();
// iterate over all possibilites for i and j, and use the formula for Matrix
// multiplication to calculate element (i,j) of the product
Matrix<T> c = Matrix<T>(rows,b.cols);
for (int i=1; i<= c.rows; ++i)
{
for (int j=1; j<=c.cols; ++j)
{
acc=T()-T();
for (int k=1; k<=cols; ++k)
{
// formula from:
// Eric W. Weisstein. "Matrix Multiplication." From MathWorld--A Wolfram
// Web Resource. [url]http://mathworld.wolfram.com/MatrixMultiplication.html[/url]
acc=acc+(*elem(i,k)*(*b.elem(k,j)));
}
*c.elem(i,j)=acc;
}
}
return c;
}
template <typename T>
Matrix<T> Matrix<T>::operator*(const T& t) const // scalar multiplication
{
Matrix<T> M(*this);
for (int r=1; r<=rows; ++r)
{
for (int c=1; c<=cols; ++c)
{
*M.elem(r,c) = (*M.elem(r,c))*t;
}
}
return M;
}
template <typename T>
Matrix<T> Matrix<T>::operator/(const T& t) const // scalar multiplication
{
Matrix<T> M(*this);
for (int r=1; r<=rows; ++r)
{
for (int c=1; c<=cols; ++c)
{
*M.elem(r,c) = (*M.elem(r,c))/t;
}
}
return M;
}
template <typename T>
bool Matrix<T>::operator==(const Matrix<T>& m) // test equality
{
if (rows!=m.numRows() || cols!=m.numCols()) {return false;} // if the
// dimensions don't match, there is no need to continue
else for (int i=1; i<=rows; ++i) // if we have gotten this far, the dimensions
// match, so now we have to check each element of both matrices
{
for (int j=1; j<=cols; ++j)
{
if (data[i][j]!=m.data[i][j]) {return false;} // If we find a position
// where THIS and m are not equal, then return false
}
}
return true; // if we have gotten this far, we have iterated through all
// elements of THIS and m without finding any inequality, so they are equal
}
template <typename T>
T* Matrix<T>::operator[](const int& index) // return a pointer to an element
// indexed in row major order
{
// change the index to (row, column) coordinates, and then call the elem
// function
int r=ceil((float)index/cols);
int c=index%cols;
if (c==0)
{
c=c+cols; // % does not work exactly like the actual modulus operator from
// number theory on negative numbers.
}
return elem(r,c);
}
template <typename T>
int Matrix<T>::length() {return rows*cols;} // The total number of elements
template <typename T>
Matrix<T> Matrix<T>::transpose() // return the transpose of THIS Matrix (swap rows and
// columns)
{
Matrix<T> m(cols,rows); // we will tempoarily store the output in m
for (int i=1; i<=rows; ++i)
{
m.setCol(i,getRow(i)); // set the ith column of m equal to the ith row
// of THIS Matrix
}
return m;
}
template <typename T>
T Matrix<T>::determinant()
{
assert(rows==cols&&rows>=2);
if (rows==2) return ((*elem(1,1))*(*elem(2,2)))-((*elem(1,2))*(*elem(2,1)));
T acc(*elem(1,1)-*elem(1,1));
for (int i=1; i<=cols; ++i)
{
Matrix<T> temp(rows-1,cols-1);
int tempR=1;
int tempC=1;
for (int r=2; r<=rows; ++r)
{
for (int c=1; c<=cols; ++c)
{
if (c!=i)
{
*temp.elem(tempR,tempC)=*elem(r,c);
++tempC;
}
}
tempC=1;
++tempR;
}
if (i%2==0) acc=acc-((*elem(1,i))*temp.determinant());
else acc=acc+((*elem(1,i))*temp.determinant());
}
return acc;
}
template <typename T>
void Matrix<T>::push_row_back(T t) // add a new row
{
data.push_back(vector<T>(cols,t)); // create a vector of all zero's and
// add it to the bottom of THIS Matrix
++rows; // increment the number if rows
}
template <typename T>
void Matrix<T>::push_col_back(T t)
{
for (int i=0; i<rows; ++i) // iterate over all rows
{
data[i].push_back(t); // add another element to the end of each row
}
++cols; // increment the number of columns
}
template <typename T>
void Matrix<T>::pop_row_back() // add a new row
{
data.pop_back(); // delete last row a row
rows--; // decrement the number if rows
}
template <typename T>
void Matrix<T>::pop_col_back() // delete the last column
{
for (int i=0; i<rows; ++i) // iterate over all rows
{
data[i].pop_back(); // delete the last element of each row
}
cols--; // decrement the number of columns
}
template <typename T>
T Matrix<T>::cofactor(const int& r, const int& c)
{
Matrix<T> M(rows-1,cols-1);
int src_row = 1;
int src_column=1;
for (int dest_row = 1; dest_row<rows; ++dest_row)
{
if (src_row==r) ++src_row; // scip row r
for (int dest_column = 1; dest_column < cols; ++dest_column)
{
if (src_column==c) ++src_column; // skip column c
*M.elem(dest_row,dest_column) = *elem(src_row,src_column);
++src_column;
}
++src_row;
src_column=1;
}
// determine sign
if ((r+c)%2 == 0) return M.determinant();
else return -M.determinant();
}
template <typename T>
Matrix<T> Matrix<T>::cofactor_matrix()
{
Matrix<T> M(rows,cols);
// iterate over all elements and compute each cofactor
for (int r = 1; r<=rows; ++r)
{
for (int c = 1; c<=cols; ++c)
{
*M.elem(r,c) = cofactor(r,c);
}
}
return M;
}
template <typename T>
inline Matrix<T> Matrix<T>::adjoint_matrix()
{
return cofactor_matrix().transpose();
}
template <typename T>
inline Matrix<T> Matrix<T>::inverse()
{
return adjoint_matrix()/determinant();
}
template <typename T>
void Matrix<T>::print() // used to print a matrix
{
cout<<endl; // start a new line
// nested for loops iterate over the entire Matrix
for (int i=1; i<=numRows(); ++i)
{
for (int j=1; j<=numCols(); ++j)
{
cout << *elem(i,j) << " "; // add element (i,j) to the ostream
}
cout<<endl; // add a new lone
}
}
#endif
|
01-Apr-2008, 14:48
2nd Week of IA Training by crystalattice
Similar Threads
Thanks
