For example, a matrix multiply algorithm depends on how the matrix is represented

(This problem will not disappear, we have to develop techniques to solve this problem)

In particular, we would like to write the Jacobi algorithm in a way that is "independent" (or appears to be independent) to how the matrix and vector are stored:

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); }

Vector sol, b; DenseMatrix P1, Q1; SparseMatrix P2, Q2; sol = Jacobi(P1, Q1, b); // Will call the "dense" // version of matrix-vector // multiplication sol = Jacobi(P2, Q2, b); // Will call the "sparse" // version of matrix-vector // multiplication

class myOldClass ----> class myNewClass:myOldClass { { virtual function1() virual function1() } } myOldClass X; myNewClass Y; \ / \ / \ / \ / \ / \ / myOldClass *ptr ptr = &X; ptr = &Y; ptr->function1(); ptr->function1();

class DenseMatrix --------------------> { public: float A[3][3]; virtual 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; } };

class SparseMatrix:DenseMatrix { public: int nnz; float val[6]; int row[6]; int col[6]; virtual 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; } };

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 = 1.0; while ( d > 0.0001 ) { cout << v1.x[0] << "\t" << v1.x[1] << "\t" << v1.x[1] << "\n" ; v2 = (*P)*v1 + (*Q)*b; d = distance(v1, v2); v1 = v2; } return(v2); }

int main(int argc, char *argv[]) { Vector sol, b; DenseMatrix P1, Q1; sol = Jacobi(&P1, &Q1, b); // Solves dense system SparseMatrix P2, Q2; sol = Jacobi(&P2, &Q2, b); // Solves sparse system }

• NOTE: function polymorphism in C++ only works when we use pointer variables.

  We must either (1) pass the parameter by reference or (2) use the C style pass address by value.

The above program uses the C style pass address by value method

Same program, but uses C++ style pass by reference (nicer)

• Extraneous variables:
  The float A[3][3]; array variable is inherited by the sparse matrix class but this variable is never used...

• Unrepresentative naming:
  The parameter type used in the function Jacobi is Jacobi(DenseMatrix *P, DenseMatrix *Q, Vector b), it creates the misunderstanding that this function will only solve dense systems... (you want to make the function as obvious as possible to minimize confusion in the future).

class Matrix { public: virtual Vector operator*(Vector &v); };
class DenseMatrix:Matrix { public: float A[3][3]; virtual 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; } };	class SparseMatrix:Matrix { public: int nnz; float val[6]; int row[6]; int col[6]; virtual 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; } };

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 = 1.0; while ( d > 0.0001 ) { cout << v1.x[0] << "\t" << v1.x[1] << "\t" << v1.x[1] << "\n" ; v2 = (*P)*v1 + (*Q)*b; d = distance(v1, v2); v1 = v2; } return(v2); }

int main(int argc, char *argv[]) { Vector sol, b; DenseMatrix P1, Q1; sol = Jacobi(&P1, &Q1, b); // Solves dense system SparseMatrix P2, Q2; sol = Jacobi(&P2, &Q2, b); // Solves sparse system }

• Note: a solution using pass-by-reference will also work. You can modify the above program and try it out yourself...