c=c(3,2,1,2,3,2,1,1,3,39,34,26)

A=matrix(c,3,4)
A

EGauss<-function(A){
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]=m
}
}
A
}

B=EGauss(A)
B

TSup<-function(B){
n=length(B[,1])
m=length(B[1,])

if ((n+1)==m){
b=B[,m]
x=0*b
x[n]= b[n]/B[n,n]
for (j in 1:(n-1)){
x[n-j]= b[n-j]
for ( i in j:1){
x[n-j]= x[n-j]-B[n-j,n-j+i]*x[n-j+i]
}
x[n-j]= x[n-j]/B[n-j,n-j]
}}

else{
b=B[,(n+1):m]
x=0*b
x[n,]= b[n,]/B[n,n]

for (j in 1:(n-1)){
x[n-j,]= b[n-j,]
for ( i in j:1){
x[n-j,]= x[n-j,]-B[n-j,n-j+i]*x[n-j+i,]
}
x[n-j,]= x[n-j,]/B[n-j,n-j]
}}
x
}

x=TSup(B)
x

A[1:3,1:3]%*%x

# Ex. 4.8.

c=c(1, 1, 1, 0, 2, 0, 1, 2, 1, 1, 2, 1, 1, 
2, 0, 3, 2, 1, 2, 1, 2, 1, 2, 3, 1, 27,23,31,31,22)

A=matrix(c,5,6)
A

B=EGauss(A)
B

x=TSup(B)
x

A[1:5,1:5]%*%x



#------------------------
c=c(3,2,1,2,3,2,1,1,3)
A=matrix(0,3,6)
A[1:3,1:3]=matrix(c,3,3)
A


A[1:3,4:6]=diag(1,3,3)

A

EGauss<-function(A){
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]=m
}
}
A
}

B=EGauss(A)

TSup<-function(B){
n=length(B[,1])
m=length(B[1,])

if ((n+1)==m){
b=B[,m]
x=0*b
x[n]= b[n]/B[n,n]
for (j in 1:(n-1)){
x[n-j]= b[n-j]
for ( i in j:1){
x[n-j]= x[n-j]-B[n-j,n-j+i]*x[n-j+i]
}
x[n-j]= x[n-j]/B[n-j,n-j]
}}

else{
b=B[,(n+1):m]
x=0*b
x[n,]= b[n,]/B[n,n]

for (j in 1:(n-1)){
x[n-j,]= b[n-j,]
for ( i in j:1){
x[n-j,]= x[n-j,]-B[n-j,n-j+i]*x[n-j+i,]
}
x[n-j,]= x[n-j,]/B[n-j,n-j]
}}
x
}

x=TSup(B)
x

x%*%A[1:3,1:3]

# Ex. 4.8.

c=c(1, 1, 1, 0, 2, 0, 1, 2, 1, 1, 2, 1, 1, 
2, 0, 3, 2, 1, 2, 1, 2, 1, 2, 3, 1, 27,23,31,31,22)

A=matrix(c,5,6)
A

B=EGauss(A)
B

x=TSup(B)
x

A[1:5,1:5]%*%x

##------------------------
#
# algoritmo para Eliminacao de Gauss com pivotamento parcial 

GaussElimPP<-function(a,b)
{
n=length(b)

A=matrix(a,n,n)
L=0*A
troca=1:n
for (i in 1:n){L[i,i]=1}

c=b

# iteracaoes para o escalonamento

for ( k in 1:(n-1) )
    {
#------------- teste pivo
c=0*c

for (i in k:n)
{
c[i]=abs(A[i,k])
}

c1=max(c)

if (c1 > 0)
    {  
        i=k  
        while ( abs(A[i,k])< c1 ) # trocando linha 
          {
               i=i+1
           }
        for ( j in 1:n )   
           {
             c[j]=A[k,j]
             A[k,j]=A[i,j]
             A[i,j]=c[j] 
            }
        d=b[k]
        b[k]=b[i]
        b[i]=d
        tr=troca[k]
        troca[k]=troca[i]
        troca[i]=tr
   } else { stop(" o sitema nao possui solucao unica")}

#------------------

      for ( i in (k+1):n )
         {
            m=A[i,k]/A[k,k]
 
            A[i,k]=0
            L[i,k]=m
            
            for ( j in (k+1): n )
               {
                  A[i,j]=A[i,j]-m*A[k,j]
   
                }

            b[i]=b[i]-m*b[k] 
         }
     }
B=matrix(0,n,2*n+2)
B[1:n,1:n]=L
B[1:n,(n+1):(2*n)]=A
B[,2*n+1]=b
B[,2*n+2]=troca
p=B
p
}

#--------------------------------------------------------------------------------------------
# Resolver o sistema triangular superior

Triansup<- function(a,b)
{
n=sqrt(length(a))

# iteracaoes

x=0*(1:n)

x[n]=b[n]/a[n,n]

for ( k in (n-1):1 )
    {
      s=0
         for ( j in (k+1): n )
            {
              s=s+a[k,j]*x[j]
            }
       x[k]=(b[k]-s)/a[k,k]
     }
x
}

#------------------------------------------------------------------------------------
# Entrada de dados 

b=c(1,2,3,4,5)
b

n=length(b)
n

a=matrix(c(54,2,85,4,1,2,34,76,0,9,5,1,11,1,1,3,5,7,81,22,4,4,8,2,78),n,n)
a
# Matriz estendida escalonada do sistema

GaussElimPP(a,b)

# Separando matriz escalonada A1 e vetor b1

B=GaussElimPP(a,b)
B

L=B[1:n,1:n]
L

U=B[1:n,(n+1):(2*n)]
U

b1=B[,2*n+1]
b1

trocalinha=B[,2*n+2]
trocalinha
# Resolvendo o sistema escalonado (triangular superior)

x=Triansup(U,b1)
print(x,digits=22)

print(solve(U,b1), digits=22)



#------------------------------------------------------------------------------------
#
# algoritmo para Eliminacao de Gauss sem pivo
#

#------------------------------------------------------------------------------------

GaussElim<-function(a,b)
{
n=length(b)

A=matrix(a,n,n)
L=0*A
for (i in 1:n){L[i,i]=1}

c=b

# iteracaoes para o escalonamento

for ( k in 1:(n-1) )
    {
      for ( i in (k+1):n )
         {
            m=A[i,k]/A[k,k]
 
            A[i,k]=0
            L[i,k]=m
            
            for ( j in (k+1): n )
               {
                  A[i,j]=A[i,j]-m*A[k,j]
   
                }

            b[i]=b[i]-m*b[k] 
         }
     }
B=matrix(0,n,2*n+1)
B[1:n,1:n]=L
B[1:n,(n+1):(2*n)]=A
B[,2*n+1]=b
p=B
p
}

#--------------------------------------------------------------------------------------------
# Resolver o sistema triangular superior

Triansup<- function(a,b)
{
n=sqrt(length(a))

# iteracaoes

x=0*(1:n)

x[n]=b[n]/a[n,n]

for ( k in (n-1):1 )
    {
      s=0
         for ( j in (k+1): n )
            {
              s=s+a[k,j]*x[j]
            }
       x[k]=(b[k]-s)/a[k,k]
     }
x
}

#------------------------------------------------------------------------------------
# Entrada de dados 

b=c(1,2,3,4,5)
b

n=length(b)
n

a=matrix(c(54,2,5,4,1,2,34,6,0,9,5,1,11,1,1,3,5,7,81,22,4,4,8,2,78),n,n)
a
# Matriz estendida escalonada do sistema

GaussElim(a,b)

# Separando matriz escalonada A1 e vetor b1

B=GaussElim(a,b)
B

L=B[1:n,1:n]
L

A1=B[1:n,(n+1):(2*n)]
A1

b1=B[,2*n+1]
b1

# Resolvendo o sistema escalonado (triangular superior)

x=Triansup(A1,b1)
print(x,digits=22)


#--------------------------------------------------
# Autovalores por LU

aut<-function(a,m){
b=0*a[1,]
n=length(a[1,])
for (i in 1:m){
B=GaussElim(a,b)
L=B[1:n,1:n]
U=B[1:n,(n+1):(2*n)]
a=U%*%L
}
for(i in 1:n){b[i]=a[i,i]}
b
}

m=100
# Teste 

a=matrix(c(1,3,1,-2),2,2)
a

print(aut(a,50),digits=22)
print(eigen(a)$values,digits=22)


#--------------------------------------------------------------------------
# Fatoração Cholesky
chol<-function(a){
n=length(a[1,])
G=0*a
G[1,1]=sqrt(a[1,1])
for ( i in 2:n){
G[i,1]=a[i,1]/G[1,1]
for(j in 2:i){ if (j< i){G[i,j]=(a[i,j]-G[i,]%*%G[j,])/G[j,j]}}
G[i,i]=sqrt(a[i,i]-G[i,]%*%G[i,])
}
G
}

C=chol(a)
C
C%*%t(C)


#------------------------------------------------------------------------
# Autovalores por QR

autc<-function(a,m){
b=0*a[1,]
n=length(a[1,])
for (i in 1:m){
B=chol(t(Conj(a))%*%a)
Q=a%*%solve(t(B))
R=t(B)
a=R%*%Q
}
for(i in 1:n){b[i]=a[i,i]}
b
}

# Teste 

a=matrix(c(54,2,5,4,1,2,34,6,0,9,5,6,11,1,1,4,10,1,81,22,1,9,1,22,78),n,n)
a # não simétrica

autc(a,200)

# Teste raízes de polinômio p(x)=x^7+x^6+x^5+x^4+x^3+x^2+x+1

a=matrix(0,7,7)
for (i in 1:6){
a[i,i+1]=1
a[7,i]=-1}
a[7,7]=-1
a # não simétrica

autc(a,200)

b=eigen(a)
b$values

# Teste raízes de polinômio p(x)=x^2-3x+2

a=matrix(c(0,-2,1,3),2,2)
a # não simétrica

autc(a,200)

b=eigen(a)
b$values