GMatrix.h

Go to the documentation of this file.

#ifndef __VECMATH_GMATRIX_HPP
#define __VECMATH_GMATRIX_HPP

#ifndef __VECMATH_MATRIX3D_HPP
#include <vecmath/Matrix3d.h>
#endif

#ifndef __VECMATH_MATRIX3F_HPP
#include <vecmath/Matrix3f.h>
#endif

#ifndef __VECMATH_MATRIX4D_HPP
#include <vecmath/Matrix4d.h>
#endif

#ifndef __VECMATH_MATRIX4F_HPP
#include <vecmath/Matrix4f.h>
#endif

#ifndef __VECMATH_GVECTOR_HPP
#include <vecmath/GVector.h>
#endif

class DoubleArray2 {
    int cols;
    bool preset;
    double *arr;
public:
    DoubleArray2() : preset(false), arr(0) {}

    DoubleArray2(int i, int j) : preset(false), arr(0) {
        resize(i, j);
    }

    ~DoubleArray2() {
        if (!preset && arr != 0) {
            //cerr << "Delete " << arr << endl;
            delete arr;
        }
    }

    void resize(int i, int j) {
        if (arr != 0)
            delete arr;
        arr = new double[i * j];
        //cerr << "New " << arr << endl;
        cols = j;
    }

    // swap arrays - used to destroy previous array on local destructor
    void swap(DoubleArray2& d) {
        Math::swap(cols, d.cols);
        Math::swap(arr, d.arr);
    }

    void setPreset(int c, double *data) {
        if (!preset && arr)
            delete arr;
        preset = true;
        cols = c;
        arr = data;
        //cout << "set preset " << arr << " cols: " << cols << endl;
    }

    double* operator [] (int i) const {
        return &(arr[i * cols]);
    }

    double& operator ()(int i, int j) const {
        return arr[i * cols + j];
    }
};


class GMatrix {
public:
        static const double EPS;// = 1E-10;
public:
    int nRow;
    int nCol;
    DoubleArray2 values;

    GMatrix(int rows, int cols) {
        values.resize(rows, cols);
        nRow = rows;
        nCol = cols;
        setIdentity();
    }

    GMatrix(int rows, int cols, const double ad[]);

    // for passing array dirrectly to this matrix - dangerous but save copying time
    GMatrix(int rows, int cols, double *data, bool preset = false);

    GMatrix(const GMatrix& gmatrix);

    virtual ~GMatrix() {}

    double & operator()(int i, int j) const {
        return values[i][j];
    }

    double * operator[](int i) const {
        return values[i];
    }


    int LUD(GMatrix& LU, GVector& permutation);
    int SVD(GMatrix& U, GMatrix& W, GMatrix& V);
    // this = U * W * transpose(V)
    void add(const GMatrix & gmatrix);
    void add(const GMatrix & gmatrix, const GMatrix & gmatrix1);
    static int computeSVD(GMatrix& gmatrix, GMatrix& gmatrix1, GMatrix& gmatrix2,
                          GMatrix& gmatrix3, bool flag);
    static int compute_2X2(double d, double d1, double d2, double ad[], double ad1[], double ad2[],
                           double ad3[], double ad4[], int i);
    static void compute_qr(int i, int j, double ad[], double ad1[], GMatrix& gmatrix,
                           GMatrix& gmatrix1, bool flag);
    static double compute_rot(double d, double d1, double ad[], double ad1[]);
    static double compute_shift(double d, double d1, double d2);

    void copySubMatrix(int i, int j, int k, int l, int i1, int j1, GMatrix& gmatrix);

    static double d_sign(double d, double d1) {
        double d2 = (d < 0.0) ? -d : d;
        return (d1 < 0.0) ? -d2 : d2;
    }

    virtual bool epsilonEquals(const GMatrix & gmatrix, double d) const;
    virtual bool equals(const GMatrix & gmatrix) const;
    void get(GMatrix& gmatrix) const;
    void get(Matrix3d& matrix3d) const;
    void get(Matrix3f& matrix3f) const;
    void get(Matrix4d& matrix4d) const;
    void get(Matrix4f& matrix4f) const;

    void getColumn(int col, GVector gvector) const {
        if (gvector.getSize() < nRow)
            gvector.setSize(nRow);
        for (int r = 0; r < nRow; r++)
            gvector.values[r] = values[r][col];
    }

    void getColumn(int col, double ad[]) const {
        for (int r = 0; r < nRow; r++)
            ad[r] = values[r][col];
    }

    double getElement(int row, int col) const {
        return values[row][col];
    }

    int getNumCol() const {
        return nCol;
    }

    int getNumRow() const {
        return nRow;
    }

    void getRow(int row, GVector gvector) const {
        if (gvector.getSize() < nCol)
            gvector.setSize(nCol);
        for (int c = 0; c < nCol; c++)
            gvector.values[c] = values[row][c];

    }

    void getRow(int row, double ad[]) const {
        for (int c = 0; c < nCol; c++)
            ad[c] = values[row][c];
    }

    virtual int hashCode() const {
        int i = 0;
//      for (int j = 0; j < nRow; j++) {
//          for (int k = 0; k < nCol; k++)
//              i ^= (new Double(values[j][k])).hashCode();
//
//      }
        return i;
    }

    void identityMinus();

    void invert() {
        invertGeneral(*this);
    }

    void invert(GMatrix& gmatrix) {
        invertGeneral(gmatrix);
    }

    static void luBacksubstitution(int i, double ad[], int ai[], double ad1[]);
    static bool luDecomposition(int i, double ad[], int ai[], int ai1[]);
    void mul(const GMatrix& gmatrix);
    void mul(const GMatrix& gmatrix, const GMatrix& gmatrix1);
    void mul(const GVector& gvector, const GVector& gvector1);
    void mulTransposeBoth(const GMatrix& gmatrix, const GMatrix& gmatrix1);
    void mulTransposeLeft(GMatrix gmatrix, GMatrix gmatrix1);
    void mulTransposeRight(const GMatrix& gmatrix, const GMatrix& gmatrix1);
    void negate();
    void negate(const GMatrix& gmatrix);
    void set(const GMatrix& gmatrix);
    void set(const Matrix3d& matrix3d);
    void set(const Matrix3f& matrix3f);
    void set(const Matrix4d& matrix4d);
    void set(const Matrix4f& matrix4f);
    void set(const double ad[]);
    void setColumn(int i, const GVector& gvector);
    void setColumn(int i, const double ad[]);

    void setElement(int i, int j, double d) {
        values[i][j] = d;
    }

    void setIdentity() {
        setScale(1.0);
    }

    void setRow(int i, const GVector& gvector);
    void setRow(int i, const double ad[]);
    void setScale(double d);
    void setSize(int i, int j);

    void setZero() {
        setZero(0, 0);
    }

    void sub(const GMatrix& gmatrix);
    void sub(const GMatrix& gmatrix, const GMatrix& gmatrix1);

//      public String toString() {
//      StringBuffer stringbuffer = new StringBuffer(nRow * nCol * 8);
//      for (int i = 0; i < nRow; i++) {
//          for (int j = 0; j < nCol; j++)
//              stringbuffer.append(values[i][j]).append(" ");
//
//          stringbuffer.append("\n");
//      } return stringbuffer.toString();
//      }

    double trace() const;
    void transpose();
    void transpose(const GMatrix& gmatrix);
//    friend ostream& operator <<(ostream& os, const GMatrix& m);

private:
    void resize(int rows, int cols) {
        if (nRow != rows || nCol != cols) {
            nRow = rows;
            nCol = cols;
            values.resize(rows, cols);
        }
    }

    void resizeSet(int rows, int cols) {
        if (nRow < rows || nCol < cols)
            resize(rows, cols);
    }

    void setZero(int r, int c) {
        for (int i = r; i < nRow; i++)
            for (int j = c; j < nCol; j++)
                values[i][j] = 0.0;
    }

    static void checkMatrix(GMatrix gmatrix);
    void invertGeneral(GMatrix& gmatrix);
    static void print_m(const GMatrix& gmatrix, const GMatrix& gmatrix1, const GMatrix& gmatrix2);
    static void print_se(double ad[], double ad1[]);
    static void print_svd(const double ad[], const double ad1[], const GMatrix& gmatrix, const GMatrix& gmatrix1);
    static void update_u(int i, GMatrix& gmatrix, double ad[], double ad1[]);
    static void update_u_split(int i, int j, GMatrix& gmatrix, double ad[], double ad1[],
                               bool flag, GMatrix& gmatrix1, GMatrix& gmatrix2);
    static void update_v(int i, GMatrix& gmatrix, double ad[], double ad1[]);
    static void update_v_split(int i, int j, GMatrix& gmatrix, double ad[], double ad1[],
                               bool flag, GMatrix& gmatrix1, GMatrix& gmatrix2);
    static void chase_across(double ad[], double ad1[], int i, GMatrix& gmatrix, bool flag);
    static void chase_up(double ad[], double ad1[], int i, GMatrix& gmatrix, bool flag);
};

#endif

---