DAPPLE Example: gauss
Source code:
// @TITLE "gauss.cc -- Gaussian elimination"
//
// Gaussian Elimination without pivoting
// written using the Data-Parallel Programming Library for Education (DAPPLE)
//
// We use parallelism in computing the multipliers, in multiplying each
// multiplier by the current row (an outer product), in subtracting
// that result from the lower right submatrix, and in updating the b matrix.
// and for doing inner products in the backsolve.
// This example makes heavy use of matrix slices, inner and outer products,
// and restricted matrix and vector contexts.
//
// David Kotz 1994
// based on gauss1.c - DFK 10/20/88
// $Id: gauss.cc,v 1.5 94/09/29 12:18:03 dfk Exp Locker: dfk $
#include
#include
#include
#include "dapple.h"
// #define DEBUG // adds some debugging output
// we'll work with NxN matrices
const int N = 4;
// which VP am I? useful for VP-specific math
const intVector VP(N, Identity);
// which row or column am I? useful for VP-specific math
const intMatrix VProw(N, N, IdentityR);
const intMatrix VPcol(N, N, IdentityC);
main()
{
// we are to solve Ax = b
floatMatrix A(N,N); // the matrix
floatVector b(N,N); // the right-hand side
floatVector x(N,N); // the variable, and our solution
floatVector correct(N,N); // the correct solution
floatVector m(N,N); // multiplier slice
int i; // loop variable
float s;
cout << "Enter " << N << "x" << N << " matrix A:" << endl;
cin >> A;
cout << "Enter " << N << " vector b:" << endl;
cin >> b;
cout << "Enter " << N << " vector x, the correct answer:" << endl;
cin >> correct;
cout << "Initial situation:" << endl;
cout << "A:" << endl;
cout << A << endl;
cout << "b:" << b << endl << endl;
// forward elimination
for (i = 0; i < N-1; i++) {
if (A[i][i] == 0) {
cout << "Can't do this matrix without pivoting" << endl;
exit(1);
}
ifp (VP > i) { // for all the rows below i
m = A[_][i] / A[i][i]; // compute multipliers using column i
// now subtract row i from all rows > i
A[_][i] = 0.0; // column i is zero, by design
}
// ... subtract row i from all rows > i
floatVector rowi(A[i][_]); // row i
ifp (VProw > i && VPcol > i) // for all rows > i, cols > i
A -= outer(m, rowi); // subtract outer product m * A[i][_]
// ... subtract row i from all rows > i
ifp (VP > i)
b -= m * b[i]; // we also do the subtraction in b
#ifdef DEBUG
cout << "After eliminating row: " << i << endl;
cout << "A:" << endl;
cout << A << endl;
cout << "b:" << b << endl;
cout << "m:" << m << endl << endl;
#endif
}
cout << "After forward elimination:" << endl;
cout << "A:" << endl;
cout << A << endl;
cout << "b:" << b << endl << endl;
// backsolve
for (i = N-1; i >= 0; i--) {
ifp (VP > i) // for columns of A, and "rows" of x, that are > i
s = b[i] - inner(A[i][_], x); // scalar minus (dot product of a row and x)
x[i] = s / A[i][i];
}
// print the answer
cout << "My answer: " << x << endl;
cout << "Correct answer: " << correct << endl;
cout << "Error vector: " << correct-x << endl;
}
Demonstration:
cat gauss.data
1 0 2 3
-1 2 2 -3
0 1 1 4
6 2 2 4
1 -1 2 1
-.185714 .228571 -.114285 .471428
gauss < gauss.data
Enter 4x4 matrix A:
Enter 4 vector b:
Enter 4 vector x, the correct answer:
Initial situation:
A:
1 0 2 3
-1 2 2 -3
0 1 1 4
6 2 2 4
b:1 -1 2 1
After forward elimination:
A:
1 0 2 3
0 2 4 0
0 0 -1 4
0 0 0 -70
b:1 0 2 -33
My answer: -0.185714 0.228571 -0.114286 0.471429
Correct answer: -0.185714 0.228571 -0.114285 0.471428
Error vector: 2.38419e-07 -4.17233e-07 7.07805e-07 -5.66244e-07