c=c(5,3,1,2,1,1,1,4,3)
A=matrix(c,3,3)
A
b=c(1,2,3);b

EGauss<-function(A,b){
n=length(A[,1])
for (j in 1:(n-1)){
for ( i in (j+1):n){
m= A[i,j]/A[j,j]
A[i,]=A[i,]-m*A[j,]
A[i,j]=0
b[i]=b[i]-m*b[j]
}
}
return(list(R=A,v= b))
}

B=EGauss(A,b);R=B$R;R;v=B$v;v


# Resolução de sistema triangular superior

TSup<-function(R,v){ # Retorna u tal que Ru=v
n=length(v);u=0*v
u[n]=v[n]/R[n,n]
for ( i in (n-1):1){
p=v[i]
for (j in (i+1):n){
p=p-u[j]*R[i,j]
}
u[i]=p/R[i,i]
}
u
}

u=TSup(R,v); u
print("Teste de erro"); A%*%u-b