GaussElim

3ff9934e-0071-4049-99d1-c9dc1e9d79f9 39 4 5 false false 0 decimal

&[DATE] &[TIME] - &[FILENAME]

GaussElim

m,n: pivot row and col r = argmax (h = m ... nr, abs(A[h, n])) (not implemented, use r=m)

M GaussElim A M : nr A rows : nc A cols : c 1 nc range : m 1 : n 1 : 16 mat m nr ≤ (n nc ≤ (& r m : A r n el 0 SameQ n n 1 + : k n 1 + (nc range : h m 1 + : h nr ≤ h h 1 + : f A h n el A m n el / : A h n el 0 : A h k el A h k el f A m k el * - : 13 mat for 21 line if m m 1 + : n n 1 + : 12 mat 31 line while A 31 line :

This for loop is because h=range(m+1,nr) fails when m=nr. Maybe can be substitute with a proper range variable.

M T RREF A M : m A rows : n A cols : k 1 : j 1 : B 0 : ro 1 m range : tmp 0 : 18 mat k m ≤ (j n ≤ (& p k : h k 1 + : h m ≤ h h 1 + : A p j el abs A h j el abs < 1 SameQ p h : continue if for T 0 ≠ (A p j el abs T < (& im k m range : A im j el 0 : j j 1 + : 13 mat jn j n range : p k ≠ tmp jn el A k jn el : A k jn el A p jn el : A p jn el tmp jn el : 13 mat 0 if f A k j el : A k jn el A k jn el f / : 12 mat B ro jn el ro k ≡ A ro jn el A ro jn el A ro j el A k jn el * - if : A ro jn el B ro jn el : k k 1 + : j j 1 + : 13 mat 51 line if 31 line while A 31 line :

Find de index of the largest element in the remainder col j The\00A0\column\00A0\is\00A0\negligible\002C\\00A0\zero\00A0\it\00A0\out Swap\00A0\k\002D\th\00A0\and\00A0\p\002D\th\00A0\rows Divide\00A0\the\00A0\pivot\00A0\row\00A0\by\00A0\the\00A0\pivot\00A0\element, and preserve A Subtract\00A0\multiples\00A0\of\00A0\the\00A0\pivot\00A0\row\00A0\from\00A0\all\00A0\the\00A0\other\00A0\rows In matlab you can use A(ro,jn)=A(ro,jn)-A(ro,j)*A(k,jn), because right side is evaluated first, but not in SMath Take the new values and go to next row and col

a b SameQ a num2str b num2str ≡ 11 line :

GaussElim Numeric

A 12345678933 mat :

A GaussElim 12303 - 6 - 00033 mat

A gausselim maple 12303 - 6 - 00033 mat

A 121 - 4 - 231 - 11 - 2 - 03 - 2234 mat :

A GaussElim 121 - 4 - 01 - 13 - 001 - 234 mat

A gausselim maple 121 - 4 - 01 - 13 - 001 - 234 mat

GaussElim Symbolic

A 1 x 345678 y 33 mat :

A GaussElim 1 x 3054 x * - 6 - 003754 x * - (* - 287 x * - (* + (* y 54 x * - (* + 54 x * - / 33 mat

A GaussElim A gausselim maple - simplify maple max 0

This seems to be a bug en maple gausselim:

A 123 x 5678 y 33 mat :

A GaussElim 123052 x * - 32 x - (* 003752 x * - (* - 62 x - (* + (* y 52 x * - (* + 52 x * - / 33 mat

A gausselim maple 12306 - 21 - y + 00323 - 8 x * + (* 5 y * + 2 x * y * - 6 / 33 mat

This is ok:

This is bad:

A GaussElim x 4 ≡ y 9 ≡ 21 sys at 12303 - 6 - 00033 mat

A gausselim maple x 4 ≡ y 9 ≡ 21 sys at 12306 - 12 - 00033 mat

RREF

A check.N A TOL RREF A rref maple - normi abs TOL < 11 line :

A check.S A 0 RREF A rref maple - simplify maple normi abs 0 ≡ 11 line :

TOL 109 -^:

Symbolic

A 12345678933 mat :

A TOL RREF 101 - 01200033 mat

A check.N 1

A 121 - 4 - 231 - 11 - 2 - 03 - 2234 mat :

A TOL RREF 1008 - 01010012 - 34 mat

A check.N 1

A 1230 x 56278 y 034 mat :

A 0 RREF 1004421 y - (358 x * - (* 642 x y * - (* - (* - 3 y 358 x * - (* 842 x y * - (* - (* + (* 358 x * - (21 y - (358 x * - (* 642 x y * - (* - (* / 0101421 y - (* 21 y - (358 x * - (* 642 x y * - (* - / 0018421 y - (358 x * - (* 642 x y * - (* - / - 34 mat

A check.S 1

Same code, but only for numerics

This is a rref function for numerical calculus only, but with only one (sometimes some eval's make the routines slower). Sure some other eval can make this faster.

M T RREF A M eval : nr A rows : nc A cols : m 1 : n 1 : B 0 : ro 1 nr range : tmp 0 : 18 mat m nr ≤ (n nc ≤ (& p m : h m 1 + : h nr ≤ h h 1 + : A p n el abs A h n el abs < p h : continue if for A p n el abs T < im m nr range : A im n el 0 : n n 1 + : 13 mat jn n nc range : p m ≠ tmp jn el A m jn el : A m jn el A p jn el : A p jn el tmp jn el : 13 mat 0 if f A m n el : A m jn el A m jn el f / : 12 mat B ro jn el ro m ≡ A ro jn el A ro jn el A ro n el A m jn el * - if : A ro jn el B ro jn el T NearZero : m m 1 + : n n 1 + : 13 mat 51 line if 31 line while A 31 line :

M rref A TOL RREF 11 line :

x T NearZero x abs T < 0 x if 11 line :

Can use this for eliminate small values with roundoff errors, or can try to not, to see where those errors appear (and make the routine faster)

A 12345678933 mat :

A TOL RREF 101 - 01200033 mat

A 121 - 4 - 231 - 11 - 2 - 03 - 2234 mat :

A TOL RREF 1008 - 01010012 - 34 mat

Time test

randi function as in Octave_R_m_n randi d_R 2 el_R 1 el - 1 + : r 1_m range : c 1_n range : a_m_n matrix : 14 mat a r c el d random 0 round_R 1 el + : 21 line :

lb 100 - :

ub 100 :

m 50 :

n 60 :

A 100 - 10012 mat m n randi : 4411 mat

A min A max 12 mat 100 - 10012 mat

t.o 0 time :

M A rref :

M 100000000000000000000000000000000000000000000000001 - 2.5 - 4.2 - 1.1 - 1.4 0.5 3.1 0.3 0.3 1.3 01000000000000000000000000000000000000000000000000 1.7 3.1 6.9 1.8 2.3 - 0.7 5.5 - 0.7 - 1.2 - 2.6 - 00100000000000000000000000000000000000000000000000 2.1 3.7 6.2 1.6 3 - 0.9 5.7 - 0.5 - 1.5 - 1.3 - 00010000000000000000000000000000000000000000000000 0.2 - 1 1.7 0 0.8 - 01 - 0.1 0.3 - 0.4 - 00001000000000000000000000000000000000000000000000 0.8 1.9 2.7 0.3 1.2 - 0.6 2.7 - 0.1 - 0.7 - 0.6 - 00000100000000000000000000000000000000000000000000 1.3 0.9 2.1 1 0.3 - 0.6 - 1.6 - 1 - 0.4 - 1.8 - 00000010000000000000000000000000000000000000000000 1.3 - 2.5 - 5.3 - 1.4 - 1.9 0.6 - 4.5 0.3 0.8 1.3 00000001000000000000000000000000000000000000000000 1.3 - 2.4 - 3.8 - 1 - 1.4 0.3 3.3 0.3 0.1 - 1.8 00000000100000000000000000000000000000000000000000 1.4 2.6 3.7 1.1 1.8 - 0.3 - 3.4 - 0.1 0.6 - 1.8 - 00000000010000000000000000000000000000000000000000 0.4 - 0.7 - 1.8 - 0.2 - 0.8 1.5 - 1.8 0.3 - 1.2 0.6 00000000001000000000000000000000000000000000000000 0.9 0.8 3.8 1.4 0.7 - 0.1 - 3.2 - 1.3 - 0.5 - 2.1 - 00000000000100000000000000000000000000000000000000 2.7 - 4.6 - 10.2 - 3 - 3.6 2.1 - 8.7 1.4 2.7 3.9 00000000000010000000000000000000000000000000000000 3.9 4.6 8.7 2 3.9 - 0.8 7.8 - 0.6 - 2.7 - 3.2 - 00000000000001000000000000000000000000000000000000 2.1 - 1.9 - 3.1 - 0.9 - 1.3 0.3 2.5 0.3 1.3 1.4 00000000000000100000000000000000000000000000000000 0.6 - 0.7 - 1.4 - 0.2 - 0.8 0 1.2 0.2 0.6 0.9 00000000000000010000000000000000000000000000000000 2.3 4.4 8.6 2.3 2.9 - 1.4 7.5 - 1.1 - 1.3 - 2.2 - 00000000000000001000000000000000000000000000000000 5.8 - 10.4 - 18.4 - 5.2 - 7.4 0.9 - 16.4 2.1 3.1 6.3 00000000000000000100000000000000000000000000000000 2.7 - 6.2 - 11.1 - 3 - 3.7 1.7 - 9.7 1.7 2.2 3.4 00000000000000000010000000000000000000000000000000 3.8 6.6 10.1 2.6 4.9 - 0.1 - 9.7 - 0.8 - 1.8 - 3.7 - 00000000000000000001000000000000000000000000000000 0.5 1.6 1.3 0.1 1 - 0.4 1.6 - 0.6 0.5 - 0.2 00000000000000000000100000000000000000000000000000 0.5 - 1.3 - 3.3 - 0.8 - 1 0.6 - 2.7 0.6 0.9 0.7 000000000000000000000100000000000000000000000000001 - 2.4 - 5.5 - 1.4 - 1.8 1.2 - 3.8 1.2 1.7 1.7 000000000000000000000010000000000000000000000000001 - 1.6 - 3.6 - 1.1 - 1.1 0.3 3.5 0.9 0.2 1.1 00000000000000000000000100000000000000000000000000 1.1 1.5 2.5 0.3 0.9 - 0.5 - 1.9 - 0.2 - 0.6 0.7 - 00000000000000000000000010000000000000000000000000 1.7 3.2 6.3 2.1 2.4 - 1.4 6.3 - 1.3 - 1.3 - 1.9 - 00000000000000000000000001000000000000000000000000 1.6 - 3.6 - 6 - 1.6 - 2.7 1.1 - 5.1 0.2 1.2 1.7 000000000000000000000000001000000000000000000000003 6.3 11.2 3.2 4.5 - 1.6 10.6 - 1.4 - 2.8 - 3.2 - 00000000000000000000000000010000000000000000000000 3.3 - 4.9 - 8.2 - 2 - 3.5 1.2 - 7.7 0.8 2.1 200000000000000000000000000001000000000000000000000 0.1 - 0.4 - 0.8 - 0.5 - 0.2 - 0.3 0.3 0.1 0.8 - 0.5 00000000000000000000000000000100000000000000000000 3.3 - 5.2 - 10.3 - 2.2 - 4.3 1.9 - 9.3 1.3 3.3 300000000000000000000000000000010000000000000000000 3.1 - 6.6 - 11.7 - 3.3 - 4.2 0.8 - 10.1 1.6 2.1 3.5 00000000000000000000000000000001000000000000000000 0.5 - 1.7 - 3.6 - 1 - 0.9 1.5 - 2.4 0.8 1.2 1.3 00000000000000000000000000000000100000000000000000 0.8 - 1.5 - 3.5 - 0.7 - 1.2 1.2 - 3.2 0.8 1.9 0.6 00000000000000000000000000000000010000000000000000 0.5 0.1 0.4 - 0.1 - 0.1 - 0 0.1 0 0.2 - 0.3 00000000000000000000000000000000001000000000000000 0.7 - 1.2 1.8 1.1 0.4 - 1 2.4 - 0.6 - 0.5 000000000000000000000000000000000000100000000000000 4.2 9.2 15.3 3.5 6.2 - 0.9 13.6 - 1.1 - 2.6 - 4.6 - 00000000000000000000000000000000000010000000000000 1.3 - 2.8 - 5.3 - 1.7 - 1.7 0.3 - 4.3 1.6 0.6 1.8 00000000000000000000000000000000000001000000000000 2.1 3 6.2 1.2 2.5 - 0.3 5.3 - 0.6 - 1.6 - 2.5 - 00000000000000000000000000000000000000100000000000 0.8 - 1.3 - 2.3 - 0.5 - 0.7 0.5 1.8 0.3 0.1 0.8 00000000000000000000000000000000000000010000000000 1.2 - 3.1 - 4.2 - 1 - 1.9 1.5 3.5 0.1 0.6 - 1.3 00000000000000000000000000000000000000001000000000 1.6 2.2 4.2 0.8 2 - 1.4 3.9 - 0.1 - 1.4 - 0.5 - 000000000000000000000000000000000000000001000000000 1.1 1.9 0.8 0.9 - 0.5 2.3 - 0.4 - 0.7 - 0.3 - 00000000000000000000000000000000000000000010000000 0.2 0.2 0.7 0.2 0.2 - 0.4 - 0.3 - 0.3 - 0.4 - 0.3 - 00000000000000000000000000000000000000000001000000 1.5 4.1 8.6 3.2 2 - 0.1 7.6 - 1.7 - 0 2.7 - 00000000000000000000000000000000000000000000100000 0.7 - 2.9 - 6.3 - 1.3 - 2.6 1.7 - 5.2 0.9 2.3 1.8 00000000000000000000000000000000000000000000010000 0.4 - 0.9 1.7 1.3 0.1 1.1 1.5 - 1.2 - 0.1 0.3 00000000000000000000000000000000000000000000001000 3.2 - 6.3 - 11.8 - 3.3 - 4.2 1.2 - 10.3 2.4 1.8 3.6 00000000000000000000000000000000000000000000000100 0.5 - 2.4 - 3.6 - 1.1 - 0.9 0.7 - 2.7 1.1 0.5 0.3 00000000000000000000000000000000000000000000000010 2.4 3.2 6.9 1.8 2.5 - 1.5 6 - 1.6 - 3.1 - 2.3 - 00000000000000000000000000000000000000000000000001 3.4 6.7 9.6 2 4.9 - 0.2 - 8.5 - 0.4 - 1.7 - 2.8 - 5060 mat

0 time t.o - 42.741 s *

A rref maple M - normi #

Enable this for check the result, but it's very slow.

Alvaro