Gauss Jordan Elimination Using DoubleArray
Last updated at 2:40 am UTC on 20 May 2009
Gauss Jordan Elimination Plugin uses the DoubleArray plugin http://wiki.squeak.org/squeak/6135. This plugin will solve a variable sized square matrix and returns the solution xVector.

GaussJordan Class

1. On the Browser Create the DoubleArray Class as a subclass of ByteArray
```        ArrayedCollection variableSubclass: #GaussJordan
instanceVariableNames: 'nRows nColumns contents'
classVariableNames: ''
poolDictionaries: ''
category: 'Gauss Jordan'```
2. We start of by writing the class side method.
```new: totalRows by: totalColumns
| instance |
instance := self new: totalRows * totalColumns.
instance setTotalRows: totalRows setTotalColumns: totalColumns.
^ instance```
3. In the instance side, a method is created to initialize the rows and columns of the matrix.
```setTotalRows: totalRows setTotalColumns: totalColumns
nRows := totalRows.
nColumns := totalColumns.
contents := self```
4. A couple of methods are created to retrieve the nRows and nColumns.
```getRows
^ nRows

getColumns
^ nColumns```
5. This following method loads the coefficients and solution bVector into an augmented matrix.
```loadAugmentedMatrix: coefficients and: bVector
"start index from 1 for easier implementation"
0
to: nRows - 1
do: [:i | 1
to: nColumns
do: [:j | j >= nColumns
ifTrue: [self
at: i * nColumns + j
put: (bVector at: i + 1)]
ifFalse: [self
at: i * nColumns + j
put: (coefficients at: i * nRows + j)]]]```
6. A method to retrieve the xVector along with a check on singular matrix.
```xVector
1
to: nRows
at: i
put: (self at: i * (nRows + 1))].
"singular matrix check"
0
to: nRows - 1
do: [:i |
singular := true.
1
to: nColumns - 1
do: [:j | (self at: i * nColumns + j)
~= 0
ifTrue: [singular := false]].
singular
ifTrue: [^ 'This is a Singular Matrix !!']].
7. Next, the main coordinator for the Gauss Jordan Elmination method is created.
```Gauss: doubleArrayTransform size: matrixRow by: matrixColumn
"Perform a Gaussian Elimination operation on a square matrix with the
result vector"
"doubleArrayTransform is only used in primitive."
| pivotRowIndex m1 |
<primitive: 'primitiveGaussMatrix' module: 'GaussFlexiPlugin'>
self halt.
m1 := doubleArrayTransform asGaussJordan: matrixRow by: matrixColumn.
1
to: nRows
do: [:i |
pivotRowIndex := m1 findPivotInColumn: i.
pivotRowIndex > i
ifTrue: [m1 swapRow: i withRow: pivotRowIndex].
m1 divideRow: i byColumn: i.
1
to: nRows
do: [:j | j = i
ifFalse: [m1 reduceRow: j byRow: i]]].
^ m1```
8. We proceed by adding the modules to complete the Gauss Jordan Class, firstly, a method to find the pivot row for each column.
```findPivotInColumn: column
"Find the row which contains the matimum value in a column and make
it a pivot row"
| pivot temp max |
max := (self at: column - 1 * nColumns + column) abs.
pivot := column.
1
to: nRows
do: [:i | (temp := (self at: i - 1 * nColumns + column) abs) > max
ifTrue: [max := temp.
pivot := i]].
^ pivot```
9. Swap the pivot row.
```swapRow: firstRow withRow: anotherRow
| a b |
a := firstRow - 1 * nColumns + 1.
b := anotherRow - 1 * nColumns + 1.
nColumns
timesRepeat: [self swap: a with: b.
a := a + 1.
b := b + 1]```
10. Divide row against the pivot row.
```divideRow: rowIndex byColumn: pivotColumn
"Pivot row by dividing own self with the pivot element of a column"
| temp a |
a := rowIndex - 1 * nColumns.
temp := self at: a + pivotColumn.
1
to: nColumns
do: [:i | temp = 0
ifFalse: [self at: a + i put: (self at: a + i)
/ temp]]```
11. Finally, reducing the pivot element in the other rows to zero.
```reduceRow: firstRowIndex byRow: pivotRowIndex
"Subtracting firstRow with the pivot row."
| a b temp |
a := firstRowIndex - 1 * nColumns.
b := pivotRowIndex - 1 * nColumns.
temp := self at: a + pivotRowIndex.
1
to: nColumns
do: [:i | self at: a + i put: (self at: a + i)
- (temp
* (self at: b + i))]```
12. This method converts the Gauss Jordan matrix into DoubleArray.
```asDoubleArray
n := self size.
1
to: n
do: [:i | answer at: i put: (self at: i) asFloat].

Part II: Writing the GaussFlexiPlugin

1. Create a plugin module that is going to store the primitive as 'GaussFlexiPlugin'
```        InterpreterPlugin subclass: #GaussFlexiPlugin
instanceVariableNames: 'nRows nCols'
classVariableNames: ''
poolDictionaries: ''
category: 'Gauss Jordan'```
2. Declare the instance variables as global variables in C.
```declareCVarsIn: cg
cg var: #nRows type: #int.
cg var: #nCols type: #in ```
3. The primitive:
```primitiveGaussMatrix
| result m1 |
self export: true.
self inline: false.
self var: #m1 type: 'double *'.
nCols := interpreterProxy stackIntegerValue: 0.
nRows := interpreterProxy stackIntegerValue: 1.
m1 := self loadArgumentMatrix: (result := interpreterProxy stackObjectValue: 2).
interpreterProxy failed
ifTrue: [^ nil].
self matrixGauss: m1.
with: m1
size: matrixRow
by: matrixColumn."
interpreterProxy pop: 3.
interpreterProxy push: result```
4. A method to load the matrix.
```loadArgumentMatrix: matrix
self returnTypeC: 'double *'.
interpreterProxy failed
ifTrue: [^ nil].
^ self
cCoerce: (interpreterProxy firstIndexableField: matrix)
to: 'double *'```
5. And finally, the method that does all the computation.
```matrixGauss: m1
"Transform matrix into a sequence of numbers"
| pivotRowIndex temp firstRowIndex anotherRowIndex rowIndex pivotMax |
self var: #m1 type: 'double *'.
self var: #pivotRowIndex type: 'int'.
self var: #temp type: 'double'.
self var: #pivotMax type: 'double'.
self var: #firstRowIndex type: 'int'.
self var: #anotherRowIndex type: 'int'.
self var: #rowIndex type: 'int'.
self var: #pivotColumn type: 'int'.
0 to: nRows - 1
do: [:i |
"Start of findPivotColumn:"
pivotMax := (m1 at: i * nCols + i) abs.
pivotRowIndex := i.
0 to: nRows - 1
do: [:j | "nRows - 1"
(temp := (m1 at: j * nCols + j) abs) > pivotMax
ifTrue: [pivotMax := temp.
pivotRowIndex := j]].
"End of findPivotColumn:"
pivotRowIndex > i
ifTrue: ["Start of swapRow: withRow:"
firstRowIndex := i * nCols.
anotherRowIndex := pivotRowIndex * nCols.
0 to: nCols - 1
do: [:a |
temp := m1 at: firstRowIndex + a.
m1 at: firstRowIndex + a
put: (m1 at: anotherRowIndex + a).
m1 at: anotherRowIndex + a put: temp
"Start of swapRow: withRow:"]].
"Start of divideRow: byColumn:"
rowIndex := i * nCols.
temp := m1 at: rowIndex + i.
0 to: nCols - 1
do: [:a | temp = 0
ifFalse: [m1 at: rowIndex + a put: (m1 at: rowIndex + a) / temp]].
"End of divideRow: byColumn:"
0 to: nRows - 1
do: [:a | a = i
ifFalse: ["Start of reduceRow: byRow:"
firstRowIndex := a * nCols.
anotherRowIndex := i * nCols.
temp := m1 at: firstRowIndex + i.
0 to: nCols - 1
do: [:b | m1 at: firstRowIndex + b put: (m1 at: firstRowIndex + b)
- (temp	* (m1 at: anotherRowIndex + b))
"End of reduceRow: byRow:"]]]]```

Part III: Testing the Gauss Jordan Elimination Plugin

1. In order to verify the accuracy of the plugin, the xVector is multiplied into the coefficients to check if the solutions correspond to the bVector. Hence, a method is created.
```multiplyWith: xVector
xVector class ~= Array
ifTrue: [^ xVector].
0
to: nRows - 1
do: [:j |
solution := 0.
1
to: nRows
do: [:i | solution := (self at: nColumns * j + i)
* (xVector at: i) + solution].
answer at: j + 1 put: solution].
2. This is a test program to run the Gauss Jordan plugin where the user has to specify only the coefficients and the bVector, this test program returns the computed bVector.
```example3a
"GaussJordan example3a."
"Prove that the flexible Gauss Jordan is actually working !!!"
"Augmented 4x4 matrix"
"2	1	-1	8	0"
"-3	-1	2	-11	0"
"-2	1	2	-3	0"
"1	2	3	4	5"
| m2 coefficientsA bVector dm2 dm3 bVectorComputed |
coefficientsA := #(2 1 -1 8 -3 -1 2 -11 -2 1 2 -3 1 2 3 4 ).
bVector := #(0 0 0 5 ).
m2 := GaussJordan new: bVector size by: bVector size + 1.
dm2 := m2 asDoubleArray.
dm3 := m2
Gauss: dm2
size: m2 getRows
by: m2 getColumns.
dm3 class ~= GaussJordan
ifTrue: [dm3 := dm3 asGaussJordan: bVector size by: bVector size + 1].
"Original matrix coefficients * Gauss Jordan xVector"
bVectorComputed := m2 multiplyWith: dm3 xVector.
^ bVectorComputed```
3. Finally, the legend for all the variables mentioned in the example.
• m2 = original augmented matrix.
• dm2 = original augmented matrix in DoubleArray.
• dm3 = Gauss Jordan solution matrix [Diagonal | xVector].
• bVectorComputed = get back the bVector.