I've found that I like to read about the things that people actually do with their lives. For example, its nice to read about what an artist is painting, or what a writer is currently writing. Granted, these things translate to a linguistic medium much easier than say, what I do. But I still like to sometimes write up what it is that I'm up to. I do this partially because its a challenge to "verbally" describe the process of what I'm doing. I also do this because I think its cool, and not too many people out there read diaries that are into this (or understand it).
My current project in CS5431 - Parallel Algorithms - is to create a program that solves a system of equations using the Jacobi Solving Algorithm on a matrix A that has been compressed using Ellpack-Ilpact compression. I'll explain.
A system of equations follows the form of something like this:
A00x0+A01x1+A02x2=b0
A10x0+A11x1+A12x2=b1
A20x0+A21x1+A22x2=b2
This is where A is a matrix of N by N dimensions holding the coefficients of the above equations, x is a matrix of N length holding the variable x to be solved for, and b is a matrix of N length holding the equations answers.
Now, there are algebraic ways to directly solve the above systems. However, most high school math teachers ignore the fact that there are also indirect methods, i.e. "lets make an educated guess, see what we get, and try again". These methods are interesting, and if nothing else, actually do get used in certain problems.
One of the more famous indirect method is the Jacobi solving algorithm. Lets go back to matrix A and rig the deck a little bit. Lets say that for any entry A[i][i], which would be a coefficient on the diagonal of A, that value is sufficiently large enough to be greater than the summation of all values in a given column of A for that row i. Lets also say that the matrix is sparse, which means that a decent number of the entries are 0. This gives us a matrix looking something like this:
100 x0 + 0 x1 + 1 x2 = b0
0 x0 + 100 x1 + .5 x2 = b1
.4 x0 + 0x1 + 100 x2 = b2
Which means A would look like:
100 0 1
0 100 .5
.4 0 100
See? Diagonally dominant. If we have a system of equations like this, than the Jacobi algorithm will produce the answers to x in a very few short iterations. But first, assume that N (the size of the dimensions of Matrix A) are large. Very large. Like 100. And if the matrix is sparse (i.e. mostly 0) then we have a whole lot of wasted space. Lets instead compress A into two matricies, called U and V. V will contain the necessary values of A. U will contain the location of those values in A. So, the dimensions of U and V will both be N rows by MAX_NON_ZEROS columns, where MAX_NON_ZEROS is of course, the maximum number of non-zeros in any row of matrix A. So in our little example above, U and V would be:
U: 0 2 V: 100 1
1 2 100 .5
0 2 .4 100
While this actually doesn't save space in our example, it does with large values of N.
So, anywho, we can use the Jacobi method now to solve for x. First, we'll define x to be x_new and x_old, and have both be 0 initially. Additionally, since we compressed the matrix, we'll say we know that the diagonal value is constant and set that equal to DIAG. This horribly simplifies the problem for us since we don't have to figure out which is the diagonal value in the compressed matricies - although we could, but its just annoying and it simplifies the example. Then, we can use the following loops to solve for x_new:
for (i = 0; i < N; i++) {
for (j = 0; j < MAX_NON_ZEROS; j++) {
x_new[i] = (1 / DIAG) * summation(b[i] - V[i][j] * x_old[U[i][j]], where i doesn't equal j)
}
}
You continue to perform the above operation until the summation of r[i]^2, where r[i] is the summation of (b[i] - V[i][j]*x_old[U[i][j]]), is sufficiently small.
See? Guess and check. Craziness.
Now, my assignment, now that we understand the above, is to take this algorithm and make it run on a parallel machine. What we do is we take a row in the matricies U and V and assign that to a particular processor. That processor than does the Jacobi algorithm on that particular row. At each iteration it shares its data, and we continue until r[i] is sufficiently small. The trick is that we have to share data on each iteration, which is costly, not to mention the synchronization efforts that have to occur between all processors. So this isn't really the best paralleled algorithm, but it is kinda cool.
I also have to do this in both UPC and MPI, two different approaches to parallelizing algorithms. I'll leave that for another talk since that gets into the difference between shared memory architecture and message passing.
I hope you enjoyed this little discussion. I know I did.
Miss Any?
Alright, we're gonna give this a shot
January 02, 2005
Let's see how badly I failed these last year
December 31, 2004
Okay, so its trendy
December 28, 2004
Its just like that asshole, Joe fucking Lieberman. Annoying and rather pointless.
December 27, 2004
Is this a typical Christmas?
December 26, 2004