# Processo de ortogonalização de Gram-Schmdit 

pe<-function(u,v){ # Produto escalar dos vetores u e v.
n=length(u)
p=u[1]*v[1]
for ( i in 2:n){p=p+u[i]*v[i]}
p
}    

mod<-function(u){sqrt(pe(u,u))}  # módulo do vetor u

proj<-function(u,v){pe(u,v)*v}  # Projeção de u sobre v (ortogonal), com norma de v igual 1

GramS<-function(A){  # Ortogonliza as colunas de A
n=length(A[1,])

Q=0*A

Q[,1]=A[,1]/mod(A[,1])

for (i in 2:n){
p=A[,i]
for (j in 1:(i-1)){
p=p-proj(A[,i],Q[,j])
}
Q[,i]=p/mod(p)
}


Q
}
#-------Teste

# Dada a matriz A e o vetor b, resolver o sistema Au=b. Escrevendo QR=A, o que implica em Ru=v, em que v=Q^*b
A=matrix(c ( 4, 6, 7, 1, 3, 4, 8, 5, 6), 3,3,byrow=TRUE); print("A");A; b=c( 3, 7, 10); print("b");b

Q=GramS(A); print("Q"); Q;

t(Q)%*%Q      # Teste de ortogonalização

R=t(Q)%*%A; print("R"); R
v=t(Q)%*%b; print("Q^*b");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