# Linear system Au=v.

c=c(1,2,1,-1,3,2,4,4,4,4,3,4,2,0,1,5)
A=matrix(c,4,4,,byrow=TRUE); A   # Matrix A
b=c( 5, 16, 22,15); b            # Vector v

f<-function(u){A%*%u-b}          # Au=v when f(u)=0.

B=t(A)%*%A; bb=t(A)%*%b           # A*Au=A*v or Bu=A*v, with B=A*A positive definite.

n=4
u=c(1,1,1,1)
                                 # Início do método
w=B%*%u-bb                        # w0
v=w                              # d1
alpha= t(w)%*%w /(t(B%*%w) %*%w); alpha=alpha[1]     # initial step size
u=u-alpha*w                       # u1

for ( j in 1:n){
w=(B%*%u-bb)
beta=t(B%*%w) %*%v/(t(B%*%v) %*%v); beta=beta[1]
v=w-beta*v                   # Euler method
alpha= t(v) %*%w/(t(B%*%v) %*%v);alpha=alpha[1]
u=u-alpha*v                      # almost Euler method
}
print(" Approximation to u");u                               
print(" Aproximation to f(u)") ;f(u) 
print(" Iteration number"); n