This is a homework for Math 630: Linear Algebra
Textbook : Linear Algebra and Its Applications, G. Strang, Thomson Brooks ISBN: 0030105676
Code Brief
- Written in Python
- All matrix indices (row and column) are Zero-based
- line 1-92 Algorithm
- line 92-332 Tests
=======================Begin of Code=============================
001 import unittest
002
003 ######################################
004 # Three elementary transformation
005 ######################################
006 def SwapRow(aMatrix, aRow1, aRow2):
007 if aRow1!=aRow2:
008 aMatrix[aRow1],aMatrix[aRow2] = aMatrix[aRow2],aMatrix[aRow1];
009
010 def ScaleRow(aMatrix, aRow, aScalar):
011 aMatrix[aRow] = map(lambda x:x*aScalar,aMatrix[aRow]);
012
013 def CombineRow(aMatrix, aAddTo, aAddBy, aScalar):
014 aMatrix[aAddTo] = map(lambda x,y:x+y*aScalar,
015 aMatrix[aAddTo],aMatrix[aAddBy]);
016
017 def DoZeroCorrection(aValue): #not tested
018 if IsZero(aValue):
019 return 0;
020 else:
021 return aValue;
022 aMatrix[aAddTo] = map(DoZeroCorrection, aMatrix[aAddTo]);
023
024 ######################################
025 # Helpers
026 ######################################
027 def GetPivotCoeff(aPivot, aTarget):
028 return -1.0*(aTarget+0.0)/(aPivot+0.0);
029
030 EPSION = 0.000001; # 1 over million
031 def IsZero(aValue):
032 return abs(aValue)<EPSION;
033
034 def IsZeroRow(aRow):
035 return len(filter(IsZero,aRow))== len(aRow);
036
037
038 ######################################
039 # Operations
040 ######################################
041 def PivotDown(aMatrix, aLastNonZeroRow, aRow, aCol):
042 for i in range(aRow+1,aLastNonZeroRow+1):
043 if IsZeroRow(aMatrix[i]):
044 continue;
045
046 CombineRow(aMatrix,i,aRow,
047 GetPivotCoeff(aMatrix[aRow][aCol],aMatrix[i][aCol]))
048
049 def MakeColumnPivot(aMatrix, aRow, aLastNonZeroRow, aCol):
050 for c in range(aCol,len(aMatrix[0])):
051 for r in range(aRow,aLastNonZeroRow+1):
052 if not IsZero(aMatrix[r][c]):
053 SwapRow(aMatrix,aRow,r);
054 return c;
055
056 # This line will never be hit, because after doing MakeRowPivot,
057 #at least, current row has non-zero entry, thus return with the
058 # shortcut in SwapRow when Row1=Row2
059 return len(aMatrix[0])-1;
060
061 def MakeRowPivot(aMatrix, aCurrRow, aLastNonZeroRow):
062 for r in range(aLastNonZeroRow,aCurrRow,-1):
063 if IsZeroRow(aMatrix[aCurrRow]):
064 SwapRow(aMatrix,aCurrRow,r);
065 else:
066 return r;
067
068 return aCurrRow;
069
070
071
072 ######################################
073 # Gaussian Elimination Algorithm
074 ######################################
075 def DoGaussianElimination(aMatrix):
076 M = len(aMatrix);
077 N = len(aMatrix[0]);
078
079 pivotCol = 0;
080 lastRow = M-1;
081 for pivotRow in range(0,M):
082 lastRow = MakeRowPivot(aMatrix, pivotRow,lastRow)
083 if pivotRow >= lastRow: #In case that all rows down are zero
084 break;
085
086 pivotCol = MakeColumnPivot(aMatrix,pivotRow,lastRow,pivotCol);
087 if pivotCol >= (N-1): #In case that all rows down are zero
088 break;
089
090 PivotDown(aMatrix,lastRow,pivotRow,pivotCol);
091
092 pivotCol +=1;
093
094
095
096 ######################################
097 # Test fixtures
098 ######################################
099 class TestEGAlgorithm(unittest.TestCase):
100 def testDoGaussianElimination_EmptyRow(self):
101 B = [ [ 0, 0]];
102 A = [ [ 0, 0]];
103 DoGaussianElimination(A);
104 self.assertEqual(A,B);
105 B = [ [ 0, 0],
106 [ 0, 0]];
107 A = [ [ 0, 0],
108 [ 0, 0]];
109 DoGaussianElimination(A);
110 self.assertEqual(A,B);
111 B = [ [ 0, 0],
112 [ 0, 0],
113 [ 0, 0]];
114 A = [ [ 0, 0],
115 [ 0, 0],
116 [ 0, 0]];
117 DoGaussianElimination(A);
118 self.assertEqual(A,B);
119
120 B = [ [11,12],
121 [ 0,-1],
122 [ 0, 0]];
123 A = [ [11,12],
124 [22,23],
125 [ 0, 0]];
126 DoGaussianElimination(A);
127 self.assertEqual(A,B);
128 A = [ [0, 0],
129 [22,23],
130 [11,12]];
131 DoGaussianElimination(A);
132 self.assertEqual(A,B);
133 A = [ [11,12],
134 [ 0, 0],
135 [22,23]];
136 DoGaussianElimination(A);
137 self.assertEqual(A,B);
138
139 B = [ [11,12],
140 [ 0, 0],
141 [ 0, 0]];
142 A = [ [ 0, 0],
143 [ 0, 0],
144 [11,12]];
145 DoGaussianElimination(A);
146 self.assertEqual(A,B);
147 A = [ [11,12],
148 [ 0, 0],
149 [ 0, 0]];
150 DoGaussianElimination(A);
151 self.assertEqual(A,B);
152
153 def testDoGaussianElimination_Rank(self):
154 # M > N case
155 A = [ [0,0],
156 [0,1],
157 [0,0],
158 [1,0]];
159 B = [ [1,0],
160 [0,1],
161 [0,0],
162 [0,0]];
163 DoGaussianElimination(A);
164 self.assertEqual(A,B);
165
166 # N > M case
167 A = [ [0,0,0,0,0,1],
168 [0,0,0,1,0,0],
169 [0,0,0,0,0,0],
170 [1,0,0,0,0,0]];
171 B = [ [1,0,0,0,0,0],
172 [0,0,0,1,0,0],
173 [0,0,0,0,0,1],
174 [0,0,0,0,0,0]];
175 DoGaussianElimination(A);
176 self.assertEqual(A,B);
177
178 def testDoGaussianElimination_Column(self):
179 A = [ [0,0,0,0,0],
180 [0,0,0,0,0],
181 [0,0,0,0,1],
182 [0,0,1,0,0],
183 [1,0,0,0,0]];
184 B = [ [1,0,0,0,0],
185 [0,0,1,0,0],
186 [0,0,0,0,1],
187 [0,0,0,0,0],
188 [0,0,0,0,0]];
189 DoGaussianElimination(A);
190 self.assertEqual(A,B);
191
192 A = [ [0,1, 3],
193 [0,2, 4]];
194 B = [ [0,1, 3],
195 [0,0,-2]];
196 DoGaussianElimination(A);
197 self.assertEqual(A,B);
198
199 A = [ [0,12,0],
200 [0,22,0],
201 [0,32,0]];
202 B = [ [0,12,0],
203 [0, 0,0],
204 [0, 0,0]];
205 DoGaussianElimination(A);
206 self.assertEqual(A,B);
207
208 A = [[0, 0, 0],
209 [21,22,23],
210 [0, 0, 0]];
211 B = [[21,22,23],
212 [0, 0, 0],
213 [0, 0, 0]];
214 DoGaussianElimination(A);
215 self.assertEqual(A,B);
216
217 def testDoGaussianElimination_Textbook(self):
218 # Simple elimination example(Textbook p.12)
219 A = [ [ 2, 1, 1, 5],
220 [ 4,-6, 0,-2],
221 [-2, 7, 2, 9]];
222 B = [ [ 2, 1, 1, 5],
223 [ 0,-8,-2,-12],
224 [ 0, 0, 1, 2]];
225 DoGaussianElimination(A);
226 self.assertEqual(A,B);
227
228 # Roundoff error test(Textbook p.62)
229 A = [ [ 0.0001, 1.0, 1],
230 [ 1.0, 1.0, 2]];
231 B = [ [ 0.0001, 1.0, 1],
232 [ 0,-9999,-9998]];
233 DoGaussianElimination(A);
234 self.assertEqual(A,B);
235
236
237 def testIsEqualToZero(self):
238 self.assertTrue(IsZero(EPSION/2.0));
239 self.assertFalse(IsZero(EPSION));
240 self.assertFalse(IsZero(EPSION+EPSION/10.0));
241
242 def testGetPivotCoeff(self):
243 self.assertEquals(GetPivotCoeff(1,2),-2);
244 self.assertEquals(GetPivotCoeff(3,1),-1.0/3.0);
245
246 def testIsZeroRow(self):
247 r = [0];
248 self.assertTrue(IsZeroRow(r));
249 r = [0,0];
250 self.assertTrue(IsZeroRow(r));
251 r = [0,0,0];
252 self.assertTrue(IsZeroRow(r));
253
254 r = [0,1];
255 self.assertFalse(IsZeroRow(r));
256 r = [0,0,1];
257 self.assertFalse(IsZeroRow(r));
258
259 r = [1,0];
260 self.assertFalse(IsZeroRow(r));
261 r = [1,0,0];
262 self.assertFalse(IsZeroRow(r));
263 r = [1,0,1];
264 self.assertFalse(IsZeroRow(r));
265
266
267 def testSwap(self):
268 a = [[11,12],[21,22]];
269 b = [[21,22],[11,12]];
270 SwapRow(a,0,1)
271 self.assertEquals(a,b);
272
273 a = [[11,12]];
274 b = [[11,12]];
275 SwapRow(a,0,0)
276 self.assertEquals(a,b);
277
278 a = [[11,12]];
279 b = [[11,12]];
280 SwapRow(a,0,0)
281 self.assertEquals(a,b);
282
283 a = [[11],[21]];
284 b = [[21],[11]];
285 SwapRow(a,0,1)
286 self.assertEquals(a,b);
287
288 a = [[11],[21]];
289 b = [[21],[11]];
290 SwapRow(a,1,0)
291 self.assertEquals(a,b);
292
293 def testScalaRow(self):
294 a = [[11,12]];
295 b = [[110,120]];
296 ScaleRow(a,0,10.0);
297 self.assertEquals(a,b);
298
299 a = [[11],[12]];
300 b = [[110],[12]];
301 ScaleRow(a,0,10.0);
302 self.assertEquals(a,b);
303
304 a = [[11,12],[21,22]];
305 b = [[110,120],[21,22]];
306 ScaleRow(a,0,10.0);
307 self.assertEquals(a,b);
308 b = [[110,120],[210,220]];
309 ScaleRow(a,1,10.0);
310 self.assertEquals(a,b);
311
312 def testCombineRow(self):
313 a = [[11],
314 [21]];
315 b = [[2111],
316 [ 21]];
317 CombineRow(a,0,1,100);
318 self.assertEquals(a,b);
319
320 a = [[11,12],
321 [21,22]];
322 b = [[2111,2212],
323 [ 21, 22]];
324 CombineRow(a,0,1,100);
325 self.assertEquals(a,b);
326
327 a = [[11,12],
328 [21,22]];
329 b = [[ 11, 12],
330 [1121,1222]];
331 CombineRow(a,1,0,100);
332 self.assertEquals(a,b);
333
334
335 if __name__ == "__main__":
336 unittest.main();
337
==========================End of Code=============================
No comments:
Post a Comment