Matrix.h

Go to the documentation of this file.

#ifndef G3D_MATRIX_H
#define G3D_MATRIX_H

#include "G3D/g3dmath.h"
#include "G3D/Vector3.h"
#include "G3D/Vector4.h"
#include "G3D/Matrix3.h"
#include "G3D/Matrix4.h"
#include "G3D/ReferenceCount.h"

namespace G3D {

class Matrix {
public:
 typedef float T;

 static int debugNumCopyOps;

 static int debugNumAllocOps;

private:
public:

 class Impl : public ReferenceCountedObject {
 public:

 static void* operator new(size_t size) {
 return System::malloc(size);
 }

 static void operator delete(void* p) {
 System::free(p);
 }

 ~Impl();

 private:
 friend class Matrix;

 T** elt;

 T* data;

 int R;

 int C;

 int dataSize;

 void setSize(int newRows, int newCols);

 inline Impl() : elt(NULL), data(NULL), R(0), C(0), dataSize(0) {}

 Impl(const Matrix3& M);

 Impl(const Matrix4& M);

 inline Impl(int r, int c) : elt(NULL), data(NULL), R(0), C(0), dataSize(0) {
 setSize(r, c);
 }

 Impl& operator=(const Impl& m);

 inline Impl(const Impl& B) : elt(NULL), data(NULL), R(0), C(0), dataSize(0) {
 // Use the assignment operator
 *this = B;
 }

 void setZero();

 inline void set(int r, int c, T v) {
 debugAssert(r < R);
 debugAssert(c < C);
 elt[r][c] = v;
 }

 inline const T& get(int r, int c) const {
 debugAssert(r < R);
 debugAssert(c < C);
 return elt[r][c];
 }

 void mul(const Impl& B, Impl& out) const;

 void add(const Impl& B, Impl& out) const;

 void add(T B, Impl& out) const;

 void sub(const Impl& B, Impl& out) const;

 void sub(T B, Impl& out) const;

 void lsub(T B, Impl& out) const;

 void arrayMul(const Impl& B, Impl& out) const;

 void mul(T B, Impl& out) const;

 void arrayDiv(const Impl& B, Impl& out) const;

 void div(T B, Impl& out) const;

 void negate(Impl& out) const;

 void inverseViaAdjoint(Impl& out) const;

 void inverseInPlaceGaussJordan();

 void adjoint(Impl& out) const;

 void cofactor(Impl& out) const;

 T cofactor(int r, int c) const;

 void transpose(Impl& out) const;

 T determinant() const;

 T determinant(int r, int c) const;

 void arrayLog(Impl& out) const;

 void arrayExp(Impl& out) const;

 void arraySqrt(Impl& out) const;

 void arrayCos(Impl& out) const;

 void arraySin(Impl& out) const;

 void swapRows(int r0, int r1);

 void swapAndNegateCols(int c0, int c1);

 void mulRow(int r, const T& v);

 void abs(Impl& out) const;

 void withoutRowAndCol(int excludeRow, int excludeCol, Impl& out) const;

 bool anyNonZero() const;

 bool allNonZero() const;

 void setRow(int r, const T* vals);

 void setCol(int c, const T* vals);
 };
private:

 typedef shared_ptr<Impl> ImplRef;

 ImplRef impl;

 inline Matrix(ImplRef i) : impl(i) {}
 inline Matrix(Impl* i) : impl(ImplRef(i)) {}

 class SortRank {
 public:
 T value;
 int col;

 inline bool operator>(const SortRank& x) const {
 return x.value > value;
 }
 
 inline bool operator<(const SortRank& x) const {
 return x.value < value;
 }
 
 inline bool operator>=(const SortRank& x) const {
 return x.value >= value;
 }
 
 inline bool operator<=(const SortRank& x) const {
 return x.value <= value;
 }
 
 inline bool operator==(const SortRank& x) const {
 return x.value == value;
 }

 inline bool operator!=(const SortRank& x) const {
 return x.value != value;
 }
 };

 Matrix vectorPseudoInverse() const;
 Matrix partitionPseudoInverse() const;
 Matrix colPartPseudoInverse() const;
 Matrix rowPartPseudoInverse() const;

 Matrix col2PseudoInverse(const Matrix& B) const;
 Matrix col3PseudoInverse(const Matrix& B) const;
 Matrix col4PseudoInverse(const Matrix& B) const;
 Matrix row2PseudoInverse(const Matrix& B) const;
 Matrix row3PseudoInverse(const Matrix& B) const;
 Matrix row4PseudoInverse(const Matrix& B) const;

public:

 Matrix() : impl(new Impl(0, 0)) {}

 Matrix(const Matrix3& M) : impl(new Impl(M)) {}

 Matrix(const Matrix4& M) : impl(new Impl(M)) {}

 template<class S>
 static Matrix fromDiagonal(const Array<S>& d) {
 Matrix D = zero(d.length(), d.length());
 for (int i = 0; i < d.length(); ++i) {
 D.set(i, i, d[i]);
 }
 return D;
 }

 static Matrix fromDiagonal(const Matrix& d);

 Matrix(int R, int C) : impl(new Impl(R, C)) {
 impl->setZero();
 }

 static Matrix zero(int R, int C);

 static Matrix one(int R, int C);

 static Matrix identity(int N);

 static Matrix random(int R, int C);

 inline int rows() const {
 return impl->R;
 }

 inline int cols() const {
 return impl->C;
 }

 Matrix& operator*=(const T& B);

 Matrix& operator/=(const T& B);

 Matrix& operator+=(const T& B);

 Matrix& operator-=(const T& B);

 Matrix& operator*=(const Matrix& B);

 Matrix& operator+=(const Matrix& B);

 Matrix& operator-=(const Matrix& B);

 Matrix subMatrix(int r1, int r2, int c1, int c2) const;

 inline Matrix operator*(const Matrix& B) const {
 Matrix C(impl->R, B.impl->C);
 impl->mul(*B.impl, *C.impl);
 return C;
 }

 inline Matrix operator*(const T& B) const {
 Matrix C(impl->R, impl->C);
 impl->mul(B, *C.impl);
 return C;
 }

 inline Matrix operator+(const Matrix& B) const {
 Matrix C(impl->R, impl->C);
 impl->add(*B.impl, *C.impl);
 return C;
 }

 inline Matrix operator-(const Matrix& B) const {
 Matrix C(impl->R, impl->C);
 impl->sub(*B.impl, *C.impl);
 return C;
 }

 inline Matrix operator+(const T& v) const {
 Matrix C(impl->R, impl->C);
 impl->add(v, *C.impl);
 return C;
 }

 inline Matrix operator-(const T& v) const {
 Matrix C(impl->R, impl->C);
 impl->sub(v, *C.impl);
 return C;
 }


 Matrix operator>(const T& scalar) const;

 Matrix operator<(const T& scalar) const;

 Matrix operator>=(const T& scalar) const;

 Matrix operator<=(const T& scalar) const;

 Matrix operator==(const T& scalar) const;

 Matrix operator!=(const T& scalar) const;

 inline Matrix lsub(const T& B) const {
 Matrix C(impl->R, impl->C);
 impl->lsub(B, *C.impl);
 return C;
 }

 inline Matrix arrayMul(const Matrix& B) const {
 Matrix C(impl->R, impl->C);
 impl->arrayMul(*B.impl, *C.impl);
 return C;
 }

 Matrix3 toMatrix3() const;

 Matrix4 toMatrix4() const;
 
 Vector2 toVector2() const;

 Vector3 toVector3() const;

 Vector4 toVector4() const;

 void arrayMulInPlace(const Matrix& B);

 void arrayDivInPlace(const Matrix& B);

 // Declares an array unary method and its explicit-argument counterpart
# define DECLARE_METHODS_1(method)\
 inline Matrix method() const {\
 Matrix C(impl->R, impl->C);\
 impl->method(*C.impl);\
 return C;\
 }\
 void method(Matrix& out) const;


 DECLARE_METHODS_1(abs)
 DECLARE_METHODS_1(arrayLog)
 DECLARE_METHODS_1(arrayExp)
 DECLARE_METHODS_1(arraySqrt)
 DECLARE_METHODS_1(arrayCos)
 DECLARE_METHODS_1(arraySin)
 DECLARE_METHODS_1(negate)

# undef DECLARE_METHODS_1

 inline Matrix operator-() const {
 return negate();
 }

 inline Matrix inverse() const {
 Impl* A = new Impl(*impl);
 A->inverseInPlaceGaussJordan();
 return Matrix(A);
 }

 inline T determinant() const {
 return impl->determinant();
 }

 inline Matrix transpose() const {
 Impl* A = new Impl(cols(), rows());
 impl->transpose(*A);
 return Matrix(A);
 }

 void transpose(Matrix& out) const;

 inline Matrix adjoint() const {
 Impl* A = new Impl(cols(), rows());
 impl->adjoint(*A);
 return Matrix(A);
 }

 Matrix pseudoInverse(float tolerance = -1) const;

 Matrix svdPseudoInverse(float tolerance = -1) const;

 inline Matrix gaussJordanPseudoInverse() const {
 Matrix trans = transpose();
 return (trans * (*this)).inverse() * trans;
 }

 void svd(Matrix& U, Array<T>& d, Matrix& V, bool sort = true) const;

 void set(int r, int c, T v);

 void setCol(int c, const Matrix& vec);

 void setRow(int r, const Matrix& vec);

 Matrix col(int c) const;

 Matrix row(int r) const;

 T get(int r, int c) const;

 Vector2int16 size() const {
 return Vector2int16(rows(), cols());
 }

 int numElements() const {
 return rows() * cols();
 }

 void swapRows(int r0, int r1);

 void swapAndNegateCols(int c0, int c1);

 void mulRow(int r, const T& v);

 bool anyNonZero() const;

 bool allNonZero() const;

 inline bool allZero() const {
 return !anyNonZero();
 }

 inline bool anyZero() const {
 return !allNonZero();
 }

 void serialize(TextOutput& t) const;

 std::string toString(const std::string& name) const;

 std::string toString() const {
 static const std::string name = "";
 return toString(name);
 }

 double normSquared() const;

 double norm() const;

 static const char* svdCore(float** U, int rows, int cols, float* D, float** V);

};

}

inline G3D::Matrix operator-(const G3D::Matrix::T& v, const G3D::Matrix& M) {
 return M.lsub(v);
}

inline G3D::Matrix operator*(const G3D::Matrix::T& v, const G3D::Matrix& M) {
 return M * v;
}

inline G3D::Matrix operator+(const G3D::Matrix::T& v, const G3D::Matrix& M) {
 return M + v;
}

inline G3D::Matrix abs(const G3D::Matrix& M) {
 return M.abs();
}

#endif