#include <stdio.h>
#include <math.h>
#define N 10
#define A {{1,4,8,4},{4,2,3,7},{8,3,6,9},{4,7,9,2}}
multmatrix(int n, double x[n][n],double y[n][n],double result[n][n]){
int i,j,k;
double soma;
for(i=0;i<n;i++){
for(j=0;j<n;j++){
soma=0.0;
for(k=0;k<n;k++)
soma += x[i][k]*y[k][j];
result[i][j]=soma;
}
}
}
imprime(int n, double x[n][n]){
int i,j;
for(i=0; i<n; i++){
for(j=0; j<n; j++)
printf("%5.2f\t",x[i][j]);
printf("\n");
}
printf("\n\n");
}
QRdecomp(int n, double a[n][n], double q[n][n], double r[n][n]){
int i,j,k;
double u[n],prod,norm;
// ui = ai - sum_j<i (qj·ai) qj
for(i=0;i<n;i++){
for(k=0;k<n;k++)
u[k] = a[k][i];
for(j=0;j<i;j++){
//produto escalar (qj·ai)
prod=0.0;
for(k=0;k<n;k++)
prod += q[k][j]*a[k][i];
r[j][i] = prod;
for(k=0;k<n;k++)
u[k] -= prod*q[k][j];
}
//qi = ui/|ui|
norm=0.0;
for(k=0;k<n;k++)
norm += u[k]*u[k];
norm = sqrt(norm);
r[i][i] = norm;
for(k=0;k<n;k++)
q[k][i]=u[k]/norm;
}
}
main(){
int i,j,k;
double a[N][N]={{0}};
double q[N][N],r[N][N]={0};
double v[N][N]={0},temp[N][N];
double epsilon=1.0e-2;
for(k=0;k<N;k++)
v[k][k]=1.0;
double kappa=6,w=2,alpha=2*kappa;
for(i=0;i<N-1;i++){
a[i][i]=alpha;
a[i][i+1]= a[i][i-1] = -kappa;
}
i=0; a[i][i]=alpha-kappa; a[i][i+1]=-kappa;
i=N-1; a[i][i]=alpha-kappa; a[i][i-1]=-kappa;
a[N-1][0] = a[0][N-1] = -kappa;
imprime(N,a);
do{
QRdecomp(N,a,q,r);
multmatrix(N,r,q,a);
multmatrix(N,v,q,r);
for(i=0;i<N;i++){
for(k=0;k<N;k++){
v[i][k]=r[i][k];
r[i][k]=0.0;
}
}
// imprime(N,a);
// imprime(N,v);
// getchar();
}while(fabs(a[0][1])>epsilon);
for (int i = 0; i < N; ++i)
{
printf("%d %f\t", i, fabs(a[i][i]));
for (int j = 0; j < N; ++j)
{
printf("%f ", v[i][j]);
}
printf("\n");
}
}
Comments