- There is an intricate dependency between algorithm and
the data that the algorithm operates on
- Takes for example the Jacobi algorithm
for solving linear equations:
A x = b | A = M - N => (M - N) x = b <=> Mx - Nx = b <=> Mx = Nx + b <=> x = M-1Nx + M-1b <=> x = Px + QbJacobi Iteration: xn+1 = Pxn + Qb - The Jacobi iteration is
in principle only dependent on
the matrix-vector multiplication and vector-vector addition
operations.
- However, depending on HOW a matrix and a vector is represented, the algorithm will be written differently !!!
- If a matrix is dense, the most obvious way
to represent the matrix is to use a
two-dimensional array.
- In that case, the matrix multiplcation
would be implemented in this way:
class Vector { public: float x[3]; Vector operator+(Vector &v) { Vector z; int i; for (i = 0; i < 3; i = i + 1) z.x[i] = x[i] + v.x[i]; return(z); } }class DenseMatrix { public: float A[3][3]; Vector operator*(Vector &v) { Vector z; int i, j; for (i = 0; i < 3; i=i+1) z.x[i] = 0; for (i = 0; i < 3; i=i+1) for (j = 0; j < 3; j=j+1) z.x[i] = z.x[i] + A[i][j]*v.x[j]; return z; } }; - And a Jacobi Iteration algorithm can be written as:
float distance(Vector a, Vector b) { float sum = 0; int i; for (i = 0; i < 3; i = i + 1) sum = sum + fabs(a.x[i] - b.x[i]); return(sum); } Vector Jacobi(DenseMatrix & P, DenseMatrix & Q, Vector b) { Vector v1, v2; int i; float d; for (i = 0; i < 3; i = i + 1) { v1.x[i] = 0; v2.x[i] = 1; } d = 0.1; while ( d > 0.0001 ) { v2 = P*v1 + Q*b; d = distance(v1, v2); v1 = v2; } return(v2); } - Example Program:
(Jacobi with dense matrix representation) ---
click here
- Recall that Sparse matrices can be stored using
the coordinate-wise system...
click here
- Illustrative Example:
[ 11 12 14 ] [ 22 23 25 ] [ 31 33 34 ] [ 42 45 46 ] [ 55 ] [ 65 66 67 ] [ 75 77 78 ] [ 87 88 ]Representation using Coordinate-wise method: Index 0 1 2 3 4 5 6 7 8 9 10 ----------------------------------------------------------- Val 11 12 14 22 23 25 31 33 34 42 45 Row 0 0 0 1 1 1 2 2 2 3 3 Col 0 1 3 1 2 4 0 2 3 1 4 Index 11 12 13 14 15 16 17 18 19 20 21 ------------------------------------------------------------ Val 46 55 65 66 67 75 77 78 87 88 - Row 3 4 5 5 5 6 6 6 7 7 - Col 5 4 4 5 6 4 6 7 6 7 -
- The matrix-vector multiplication using the Coordinate-wise method
is implemented as follows:
Multiply matrix (stored in Coordinate-wise method) with vector d[N] for (k = 0; k < N; k = k + 1) result[i] = 0; for (k = 0; k < nnz; k = k + 1) result[Row[k]] = result[Row[k]] + Val[k]*d[Col[k]]; - In that case, the matrix multiplcation
would be implemented in this way:
class Vector { public: float x[3]; Vector operator+(Vector &v) { Vector z; int i; for (i = 0; i < 3; i = i + 1) z.x[i] = x[i] + v.x[i]; return(z); } }class SparseMatrix { public: int nnz; float val[6]; int row[6]; int col[6]; Vector operator*(Vector &v) { Vector z; int k; for (k = 0; k < 3; k = k + 1) z.x[k] = 0; for (k = 0; k < nnz; k = k + 1) z.x[row[k]] = z.x[row[k]] + val[k]*v.x[col[k]]; return z; } }; - And we
must write another Jacobi function:
Vector Jacobi(SparseMatrix & P, SparseMatrix & Q, Vector b) { Vector v1, v2; int i; float d; for (i = 0; i < 3; i = i + 1) { v1.x[i] = 0; v2.x[i] = 1; } d = 0.1; while ( d > 0.0001 ) { v2 = P*v1 + Q*b; d = distance(v1, v2); v1 = v2; } return(v2); }We CANNOT use the original Jacobi function, because the input parameters are of DIFFERENT TYPES
- Example Program:
(Jacobi with sparse matrix representation) ---
click here
- A more elegant approach is a solution where
the algorithm is independent on
the underlying representation of the data
- Clearly, it is impossible to develop algorithms that are
truly independent on the underlying representation of the data
that the algorithm manipulates.
BUT, it is possible to adapt the data manipulation algorithm automatically
- The idealized solution will look as follows:
Vector Jacobi(Matrix P, Matrix Q, Vector b) { Vector v1, v2; int i; float d; for (i = 0; i < 3; i = i + 1) { v1.x[i] = 0; v2.x[i] = 1; } d = 0.1; while ( d > 0.0001 ) { v2 = P*v1 + Q*b; d = distance(v1, v2); v1 = v2; } return(v2); }This function implements the Jacobi algorithm in general, using the matrix-vector multiplication and vector-vector addition operations.
- And we would be able to use
the SAME function (without any change in parameter types)
to solve linear equations where the matrix is represented either
with the dense or the sparse (or any other) representation:
DenseMatrix P1, Q1; Vector sol, b; sol = Jacobi(P1, Q1, b); SparseMatrix P2, Q2; sol = Jacobi(P2, Q2, b);
- This is indeed possible with object oriented programming techniques, but before we can explain the technique, we need some background knowledge on derived classes and polymorphism.