next up previous contents
Next: Accessing the matrix size Up: Design and implementation Previous: Matrix declaration

Matrix algebra

In symbolic matrix algebra, we use expressions like C = A + B, C = A*B and tex2html_wrap_inline2547 . If an algorithm could be presented to a compiler in this way, it would make the program easier to read, and easier to maintain. The numerical package MatLab [] is already doing this, and so, to some extent, is High Performance Fortran. With C++ it is possible to implement this kind of syntax for almost any structure, although one should be careful when doing so.

The mathematical operations of interest to matrix algebra are the assignment (=), addition (+), subtraction (-) and multiplication (*) operators. There are also the update operators tex2html_wrap_inline2557 , tex2html_wrap_inline2559 and tex2html_wrap_inline2561 which can be used instead of C = C + A, C = C - A and C = C * A respectively.

The division operator might have been used for solving equations, but there is the problem of stability for this operation, which depends on the values of the matrix. Such solving is further made difficult because the solving may be done from the left side ( tex2html_wrap_inline2569 ), or the right side ( tex2html_wrap_inline2571 ). Also, solving systems usually involves the selection of a method suited to the problem, and inversions of a matrix is both time consuming and may be numerically unstable. For these reasons, it is better to specify the method explicitly in the application. How SymPLA handles this is discussed in Section gif. It might also be feasible to let division be performed as an implicit QR solver step, or to let the algorithm select a method as the default solver method.

The only division operator supported by SymPLA is used to divide all elements of the matrix with a scalar value. The update division tex2html_wrap_inline2573 is also supported for this method of scaling.

The arithmetic operations on matrices do, of course, check that the operations are made between operands of compatible dimensions, and, if that is not the case, report the error through the C++ exception handling system.

The transpose of a matrix, tex2html_wrap_inline2495 , is often used in expressions, and a number of libraries have used a syntax like A.t(), A.trans() or transp(A) to get the transpose of the matrix A. While this kind of function-call may or may not involve duplication of data, this syntax gives the impression that to use the transpose of a matrix you have to use a function with duplication of data first. In LAPACK++, the use of implicit transpose is supported by functions like Blas_Mat_Trans_Mat_Mult for multiplication involving one transposed matrix ( tex2html_wrap_inline2579 ), while in BLAS this is specified with flags.

In SymPLA, the transpose of a matrix is considered to be an attribute of the matrix, that may be accessed by A.T, which might (theoretically) have been implemented with duplication of data by a conversion operator, but instead the same class is handling transposed and nontransposed matrices. A flag in the class is used to decide whether the operation shall use the transposed or nontransposed matrix in its execution, and optimization of the operation is left to the internal routines in the library. Using this convention tex2html_wrap_inline2579 can be written A * B.T.

  enum MatrixTranspose {NONTRANSPOSED, TRANSPOSED};

  class Matrix;

  class MatrixHandler{
    private:
      MatrixTranspose transposed;

    public:
      Matrix operator +(MatrixHandler &);
      Matrix operator /(double);

      MatrixHandler &operator =(MatrixHandler &);
      MatrixHandler &operator +=(MatrixHandler &);
      MatrixHandler &operator /=(double);
  };

  class Matrix : public MatrixHandler{
    public:
      MatrixHandler T;

      Matrix &operator =(Matrix &);
      Matrix &operator =(MatrixHandler &);
  };

The skeleton class definition above shows how operations and handling of transposed matrices are implemented in this library.


next up previous contents
Next: Accessing the matrix size Up: Design and implementation Previous: Matrix declaration