Squeak
  links to this page:    
View this PageEdit this PageUploads to this PageHistory of this PageTop of the SwikiRecent ChangesSearch the SwikiHelp Guide
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.


Part I: Writing the GaussJordan Class


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
    	| answer singular |
    	answer := Array new: nRows.
    	1
    		to: nRows
    		do: [:i | answer
    				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 !!']].
    	^ answer
  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 answer |
    	n := self size.
    	answer := DoubleArray new: n.
    	1
    		to: n
    		do: [:i | answer at: i put: (self at: i) asFloat].
    	^ answer


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.
    	"matrixAddComposeMatrix: m1
    	with: m1
    	size: matrixRow
    	by: matrixColumn."
    	interpreterProxy pop: 3.
    	interpreterProxy push: result
  4. A method to load the matrix.
    loadArgumentMatrix: matrix 
    	"Load the argument 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 
    	| answer solution |
    	xVector class ~= Array
    		ifTrue: [^ xVector].
    	answer := Array new: nRows.
    	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].
    	^ answer
  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.
    	m2 loadAugmentedMatrix: coefficientsA and: bVector.
    	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.