> restart;
> with(LinearAlgebra):
> # Vorwaertselimination, ueberschreibt b;
> # Annahme: Diagonale 1 !!;
> vorelim1:=proc(L,b) 
> local n,j,k; 
> n:=Dimension(b); 
> for j from 1 by 1 to n do 
>    for k from 1 by 1 to j-1 do b[j]:=b[j]-L[j,k]*b[k] end do
> end do;
> end proc:
> #Rueckwaertselimination, ueberschreibt b;
> ruecksubs:=proc(U,b) 
> local n, j, k; 
> n:=Dimension(b); 
> for j from n by -1 to 1 do 
>    for k from n by -1 to j+1 do b[j]:=b[j]-U[j,k]*b[k] end do; 
>    b[j]:=b[j]/U[j,j] 
> end do; 
> end proc:
> #Gauss-Elimination mit Ueberschreiben von A, mit Spaltenpivot;
> gaussSP:=proc(A) 
> global p; 
> local n,m,piv,pivj,j,k,x; 
> #option trace;
> n:=RowDimension(A); 
> p:=Vector(n, i->i); 
> for m from 1 by 1 to n do  
>    piv:= abs(A[m,m]); pivj:= m; 
>    for j from m+1 by 1 to n do 
>        if abs(A[j,m]) > piv then piv := abs(A[j,m]); pivj:=j end if
>    end do; 
>    k:=p[m]; p[m]:=p[pivj]; p[pivj] :=k; 
>    for k from 1 by 1 to n do x[k]:=A[m,k]; A[m,k]:=A[pivj, k]; A[pivj,k]:= x[k] end do; 
>    for j from m+1 by 1 to n do 
>        A[j,m]:=A[j,m]/A[m,m];  
>        for k from m+1 by 1 to n do A[j,k]:=A[j,k]- A[j,m]*A[m,k] end do;
>    end do;
> end do;  
> end proc:
> # Loesen eines Gleichungssystems (mit Spaltenpivot);
> loesAxb1:=proc(A,b) 
> local j, bb,n; 
> #option trace;
> n:=Dimension(b);
> bb:=Vector(n); 
> gaussSP(A); 
> for j from 1 by 1 to n do 
>     bb[j]:=b[p[j]] 
> end do; 
> vorelim1(A,bb); 
> ruecksubs(A,bb); 
> print(bb) 
> end proc;
> A:=Matrix([[6.,7.,0,1,2,5],[13,4,2,4,5,3],[6,2,3,6,8,9],[1,3,0,0,2,4],[1,1,3,5,7,9],[0,0,4,5,0,2]]);
> b:=Vector([1,0,1,0,1,0]);
> loesAxb1(A,b);
> #Kontrolle
> A:=Matrix([[6.,7.,0,1,2,5],[13,4,2,4,5,3],[6,2,3,6,8,9],[1,3,0,0,2,4],[1,1,3,5,7,9],[0,0,4,5,0,2]]);
> b:=Vector([1,0,1,0,1,0]);
> 
> x:=LinearSolve(A,b);
> 
>