[GSoC] Week 1: Implemented Basic Matrix Classes

Hi All, When the first week of the GSoC comes to an end, I was able to finish the basic structure of the Matrix class and implement some of the functionalities related to the `DenseMatrix` class. I implemented Matrix addition, multiplication, scalar addition and scalar multiplication. I am yet to update the PR with these changes. I am currently working on Gaussian elimination and having a few problems with that. I wish to solve these as soon as possible. After sorting it out, I wish to benchmark my implementation against other libraries such as Linbox and GiNaC.

Matrix Class Structure

  1. `MatrixBase` Class
    class MatrixBase: public Basic {
    public:
        MatrixBase(unsigned row, unsigned col)
            : row_{row}, col_{col} {};    // Below methods should be implemented by the derived classes. If not
        // applicable, raise an exception    // Get and set elements
        virtual RCP<const Basic>get(unsigned i) const = 0;
        virtual void set(unsigned i, RCP<const Basic> &e) = 0;

        virtual unsigned rank() const = 0;
        virtual RCP<const Basic> det() const = 0;
        virtual RCP<const MatrixBase> inv() const = 0;

        // These functions create a new instance of either DenseMatrix or
        // SparseMatrix and return a reference to the result
        virtual RCP<const MatrixBase> add_matrix(const MatrixBase &other) const = 0;
        virtual RCP<const MatrixBase> mul_matrix(const MatrixBase &other) const = 0;

    protected:
        // Stores the dimension of the Matrix
        unsigned row_;
        unsigned col_;
    };

    Most important thing about the `MatrixBase` class is that it is derived from the CSymPy `Basic` class. This enables the classes that are derived from `MatrixBase` to be used wherever a `Basic` class is used. This has a lot of advantages. `DenseMatrix` class is derived from `MatrixBase` class.

  2. `DenseMatrix` class
    class DenseMatrix: public MatrixBase {
    public:
        // Constructors
        DenseMatrix(unsigned row, unsigned col);
        DenseMatrix(unsigned row, unsigned col, std::vector<RCP<const Basic>> &l);    // Virtual functions inherited from Basic class
        virtual std::size_t __hash__() const;
        virtual bool __eq__(const Basic &o) const;
        virtual int compare(const Basic &o) const;    // Should implement all the virtual methods from MatrixBase
        // and throw an exception if a method is not applicable.
        // add_matrix, mul_matrix will have to find the correct function
        // to call depending on the `other` argument.

        // Get and set elements
        virtual RCP<const Basic> get(unsigned i) const;
        virtual void set(unsigned i, RCP<const Basic> &e);

        virtual unsigned rank() const;
        virtual RCP<const Basic> det() const;
        virtual RCP<const MatrixBase> inv() const;

        // Matrix addition
        virtual RCP<const MatrixBase> add_matrix(const MatrixBase &other) const;
        friend RCP<const DenseMatrix> add_dense_dense(const DenseMatrix &A,
                const DenseMatrix &B);
        friend RCP<const DenseMatrix> add_dense_scalar(const DenseMatrix &A,
                RCP<const Basic> &k);
        
        // Matrix multiplication
        virtual RCP<const MatrixBase> mul_matrix(const MatrixBase &other) const;
        friend RCP<const DenseMatrix> mul_dense_dense(const DenseMatrix &A,
                const DenseMatrix &B);
        friend RCP<const DenseMatrix> mul_dense_scalar(const DenseMatrix &A,
                RCP<const Basic> &k);
        
        // Gaussian elimination
        friend RCP<const DenseMatrix> gaussian_elimination(const DenseMatrix &A);
                
    protected:
        // Matrix elements are stored in row-major order
        std::vector<RCP<const Basic>> m_;
    };

    Some of the functions are declared as friend functions so that these can access private and protected members of the `DenseMatrix` class. This makes coding the algorithms a bit easier

  3. Friend function implementations
    • add_dense_dense
      RCP<const DenseMatrix> add_dense_dense(const DenseMatrix &A,
              const DenseMatrix &B)
      {
          unsigned row = A.row_;
          unsigned col = A.col_;    CSYMPY_ASSERT(row == B.row_ && col == B.col_)    std::vector<RCP<const Basic>> sum(row*col);

          std::vector<RCP<const Basic>>::const_iterator ait = A.m_.begin();
          std::vector<RCP<const Basic>>::const_iterator bit = B.m_.begin();

          for(auto &it: sum) {
              it = add(*ait, *bit);
              ait++;
              bit++;
          }

          return rcp(new DenseMatrix(row, col, sum));
      }

    • mul_dense_dense
      RCP<const DenseMatrix> mul_dense_dense(const DenseMatrix &A,
              const DenseMatrix &B)
      {
          unsigned row = A.row_;
          unsigned col = A.col_;    CSYMPY_ASSERT(col == B.row_)

          std::vector<RCP<const Basic>> prod(row*B.col_);

          for (unsigned r = 0; r<row; r++) {
              for (unsigned c = 0; c<B.col_; c++) {
                  prod[r*B.col_ + c] = zero; // Integer Zero
                  for (unsigned k = 0; k<col; k++)
                      prod[r*B.col_ + c] = add(prod[r*B.col_ + c],
                              mul(A.m_[r*col + k], B.m_[k*B.col_ + c]));
              }
          }

          return rcp(new DenseMatrix(row, B.col_, prod));
      }

Plan for this week

I wish to finish implementing Gaussian elimination and do some benchmarks.

 

 

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: