The header file used by Matrix.cpp

Matrix.h (text/plain), 9.38 KB, created by Damjan Jovanovic on 2004-02-27 02:05:33 UTC

typedef unsigned long u32;
typedef signed long s32;

// Forward declaration
template<typename T>
class Matrix;

// Expression template classes:

template<typename T, typename Left, typename Right>
class MatrixAddPair
{
	const Left   &left;
	const Right  &right;

  public:
	MatrixAddPair(const Left &l, const Right &r) : left(l), right(r) {}

	T operator()(u32 i, u32 j) const {
		return left(i,j) + right(i,j);
	}
	u32 Columns() const {return left.Columns();}
	u32 Rows() const {return left.Rows();}
};

template<typename T, typename Left, typename Right>
class MatrixSubPair
{
	const Left   &left;
	const Right  &right;

  public:
	MatrixSubPair(const Left &l, const Right &r) : left(l), right(r) {}

	T operator()(u32 i, u32 j) const {
		return left(i,j) - right(i,j);
	}
	u32 Columns() const {return left.Columns();}
	u32 Rows() const {return left.Rows();}
};

template<typename T, typename Left, typename Right>
class MatrixMulPair
{
	const Left   &left;
	const Right  &right;

  public:
	MatrixMulPair(const Left &l, const Right &r) : left(l), right(r) {}

	T operator()(u32 i, u32 j) const {
		T sum = 0;
		const u32 rows = right.Rows();
		for (u32 k = 0; k < rows; k++)
			sum += left(i,k)*right(k,j);
		return sum;
	}
	u32 Columns() const {return right.Columns();}
	u32 Rows() const {return left.Rows();}
};

template<typename T, typename M, typename Unused=int>
class MatrixMulScalarPair
{
	const M	 &matrix;
	const T   scalar;

  public:
	MatrixMulScalarPair(const M &m, T t) : matrix(m), scalar(t) {}

	T operator()(u32 i, u32 j) const {
		return matrix(i,j)*scalar;
	}
	u32 Columns() const {return matrix.Columns();}
	u32 Rows() const {return matrix.Rows();}
};


// Operators

//Addition
template<typename T,
		 typename L1, typename L2,
		 typename R1, typename R2,
		 template<typename,typename,typename> class Left,
		 template<typename,typename,typename> class Right>
MatrixAddPair<T, Left<T,L1,L2>, Right<T,R1,R2> >
inline operator+(const Left<T,L1,L2> &l, const Right<T,R1,R2> &r)
{
	return MatrixAddPair<T, Left<T,L1,L2>, Right<T,R1,R2> >   (l,r);
}

template<typename T,
		 typename L1, typename L2,
		 template<typename,typename,typename> class Left>
MatrixAddPair<T, Left<T,L1,L2>, Matrix<T> >
inline operator+(const Left<T,L1,L2> &l, const Matrix<T> &r)
{
	return MatrixAddPair<T, Left<T,L1,L2>, Matrix<T> >   (l,r);
}

template<typename T,
		 typename R1, typename R2,
		 template<typename,typename,typename> class Right>
MatrixAddPair<T, Matrix<T>, Right<T,R1,R2> >
inline operator+(const Matrix<T> &l, const Right<T,R1,R2> &r)
{
	return MatrixAddPair<T, Matrix<T>, Right<T,R1,R2> >   (l,r);
}

template<typename T>
MatrixAddPair<T, Matrix<T>, Matrix<T> >
inline operator+(const Matrix<T> &l, const Matrix<T> &r)
{
	return MatrixAddPair<T, Matrix<T>, Matrix<T> >   (l,r);
}

// Subtraction
template<typename T,
		 typename L1, typename L2,
		 typename R1, typename R2,
		 template<typename,typename,typename> class Left,
		 template<typename,typename,typename> class Right>
MatrixSubPair<T, Left<T,L1,L2>, Right<T,R1,R2> >
inline operator-(const Left<T,L1,L2> &l, const Right<T,R1,R2> &r)
{
	return MatrixSubPair<T, Left<T,L1,L2>, Right<T,R1,R2> >   (l,r);
}

template<typename T,
		 typename L1, typename L2,
		 template<typename,typename,typename> class Left>
MatrixSubPair<T, Left<T,L1,L2>, Matrix<T> >
inline operator-(const Left<T,L1,L2> &l, const Matrix<T> &r)
{
	return MatrixSubPair<T, Left<T,L1,L2>, Matrix<T> >   (l,r);
}

template<typename T,
		 typename R1, typename R2,
		 template<typename,typename,typename> class Right>
MatrixSubPair<T, Matrix<T>, Right<T,R1,R2> >
inline operator-(const Matrix<T> &l, const Right<T,R1,R2> &r)
{
	return MatrixSubPair<T, Matrix<T>, Right<T,R1,R2> >   (l,r);
}

template<typename T>
MatrixSubPair<T, Matrix<T>, Matrix<T> >
inline operator-(const Matrix<T> &l, const Matrix<T> &r)
{
	return MatrixSubPair<T, Matrix<T>, Matrix<T> >   (l,r);
}

// Multiplication
template<typename T,
		 typename L1, typename L2,
		 typename R1, typename R2,
		 template<typename,typename,typename> class Left,
		 template<typename,typename,typename> class Right>
MatrixMulPair<T, Left<T,L1,L2>, Right<T,R1,R2> >
inline operator*(const Left<T,L1,L2> &l, const Right<T,R1,R2> &r)
{
	return MatrixMulPair<T, Left<T,L1,L2>, Right<T,R1,R2> >   (l,r);
}

template<typename T,
		 typename L1, typename L2,
		 template<typename,typename,typename> class Left>
MatrixMulPair<T, Left<T,L1,L2>, Matrix<T> >
inline operator*(const Left<T,L1,L2> &l, const Matrix<T> &r)
{
	return MatrixMulPair<T, Left<T,L1,L2>, Matrix<T> >   (l,r);
}

template<typename T,
		 typename R1, typename R2,
		 template<typename,typename,typename> class Right>
MatrixMulPair<T, Matrix<T>, Right<T,R1,R2> >
inline operator*(const Matrix<T> &l, const Right<T,R1,R2> &r)
{
	return MatrixMulPair<T, Matrix<T>, Right<T,R1,R2> >   (l,r);
}

template<typename T>
MatrixMulPair<T, Matrix<T>, Matrix<T> >
inline operator*(const Matrix<T> &l, const Matrix<T> &r)
{
	return MatrixMulPair<T, Matrix<T>, Matrix<T> >   (l,r);
}

// Matrix * scalar
template<typename T,
		 typename L1, typename L2,
		 template<typename,typename,typename> class Left>
MatrixMulScalarPair<T, Left<T,L1,L2> >
inline operator*(const Left<T,L1,L2> &l, T r)
{
	return MatrixMulScalarPair<T, Left<T,L1,L2> >   (l,r);
}

template<typename T>
MatrixMulScalarPair<T, Matrix<T> >
inline operator*(const Matrix<T> &l, T r)
{
	return MatrixMulScalarPair<T, Matrix<T> >   (l,r);
}

// Scalar * matrix
template<typename T,
		 typename R1, typename R2,
		 template<typename,typename,typename> class Right>
MatrixMulScalarPair<T, Right<T,R1,R2> >
inline operator*(T r, const Right<T,R2,R2> &l)
{
	return MatrixMulScalarPair<T, Right<T,R1,R2> >   (l,r);
}

template<typename T>
MatrixMulScalarPair<T, Matrix<T> >
inline operator*(T r, const Matrix<T> &l)
{
	return MatrixMulScalarPair<T, Matrix<T> >   (l,r);
}

// The matrix class
template<typename T>
class Matrix
{
    u32  rows;
    u32  columns;
    T    **entry;
    
  public:
    // Constructors
    Matrix(const Matrix&);
    Matrix(u32 newRows, u32 newColumns);

    // Comparison and assignment
	template<typename RHS>
    Matrix&       operator=(const RHS&);
	template<typename RHS>
    bool          operator==(const RHS&) const;

	// Mutating operators on this matrix
    Matrix&       operator+=(const Matrix&);
    Matrix&       operator-=(const Matrix&);
    Matrix&       operator*=(const T);

    // Inlines
	T operator()(u32 i, u32 j) const {
		return entry[i][j];
	}
    T& operator()(u32 i, u32 j) {
        return entry[i][j];
    }
    u32    Columns() const {return columns;}
    u32    Rows() const {return rows;}

	// Static methods

	// Householder matrix for column vector x
	static Matrix Householder(Matrix x) {
		if (x.Columns() != 1) throw SizeMismatch;
		T alpha = sqrt((x.Transposed() * x)(0,0));
		if (x(0,0) >= 0) alpha = -alpha;
		x(0,0) -= alpha;
		x *= 1/sqrt((x.Transposed() * x)(0,0));
		Matrix M(x.Columns(), x.Columns());
		M.LoadIdentity();
		M = M - 2.0f*x*x.Transposed();
		return M;
	}

    // Getters:

    // Adjoint
    Matrix Adjoint() const;

    // Determinant (by Gaussian elimination)
    T      Det() const;

    // Determinant (by Laplace expansion - slow)
    T      Det2() const;

    // Displays using C++ streams
    void   Display() const;

    // Is it diagonally dominant?    
    bool   IsStrictlyDiagonallyDominant() const;

    // l1 norm
    T      L1Norm() const;

    // l infinity norm
    T      LInfNorm() const;

    // Submatrix without a row and column
    Matrix SubmatrixWithout(u32 row, u32 col) const;

    // Returns the transpose
    Matrix Transposed() const;



    // Setters

    // Divides row by scalar
    void DivideRow(u32, T);

    // Inverts using Gauss-Jordan reduction with
    // partial pivoting
    bool Invert();
    
    // Inverts using 1/det(a)*Adj(A)
    bool Invert2();
    
    // Makes matrix into I
    void LoadIdentity();
    
    // Fills with 0's
    void LoadZeroes();

    // Swaps 2 rows
    void SwapRows(u32, u32);

    // Row subtraction
    void SubRows(u32 source, T multiple, u32 dest,
                 u32 start, u32 end);

    // Tranposes the matrix
    void Transpose();

    ~Matrix();
};

class SizeMismatch{};

// Functions

// Gauss-Seidel solution of STRICTLY
// DIAGONALLY DOMINANT matrix, with error less
// than err, in at most maxIterations iterations
template<typename T>
void GaussSeidel(const Matrix<T> &A, const Matrix<T> &b,
                 Matrix<T> &x, T err,
                 u32 maxIterations = 100);

// Gives the in-place LU decomposition of A
// with rowSwap having at least as many elements as
// A's rows. Returns true on success, false otherwise.
template<typename T>
bool LUDecompose(Matrix<T> &A, u32 rowSwap[]);

// Solves for the column vector b, in place, using
// the LU decomposition of A given by LUDecompose()
template<typename T>
void LUSolve(const Matrix<T> &A, Matrix<T> &b,
             const u32 rowSwap[]);

// Positive integer power of a matrix
template<typename T>
Matrix<T> pow(const Matrix<T> &A, u32 exponent);

// Solution of Ax = b by Gaussian elimination.
// Matrix A may be any size that matches b.
// A and b are always changed.
// If unique solution, Solve() return 1,
// b is changed to x.
// If no solution, Solve() return 0.
// If infinite solutions, basis contains vectors
// of the form form
// [x1 x2 x3 x4 ... xn], so that
// x = x1 + x2*s + x3*t + ... + xn*u,
// x's column vectors; s,t...u parameters;
// n is returned.
template<typename T>
u32 Solve(Matrix<T> &A, Matrix<T> &b, Matrix<T> *basis = 0);