>> help gauss_elim x = gauss_elim(A,b,i_test) takes a square matrix A and a vector b and returns the solution of Ax=b. If i_test = 0 then there is pivoting, otherwise there is none. NOTE: the routine writes over A and b, so you can't use it twice in a row! >> A = rand(100,100); b = rand(100,1); >> A_o = A; b_o = b; No pivoting >> x = gauss_elim(A,b,1); >> max(abs(A*x-b)) ans = 5.8917e-11 >> A = A_o; b = b_o; With pivoting >> x = gauss_elim(A,b,0); >> max(abs(A*x-b)) ans = 2.5441e-13 Notice that the residual is 100 times smaller with pivoting! to see an operation count, do the following. "flops(0)" zeros out the flops counter. At the end of each stop, I find the number of flops and put that into counter. Flops only counts arithmetic operations. >> for n=1:8, flops(0); nn = 2^n; yy(n) = nn; A = rand(nn,nn); b = rand(nn,1); x = gauss_elim(A,b,1); counter(n) = flops; end note: one could fit counter to find exactly what polynomial of n the routine is taking. >> plot(yy,counter) >> hold on >> plot(yy,counter,'o') >> diary off