Fri. Dec 8th, 2023

# Calculation of the inverse of a matrix A after performing a QR factorization

``````//Calculul inversei unei matrice A dupa realizarea unei factorizari QR
#include <iostream.h>
#include <stdlib.h>
#include <conio.h>
#include <math.h>
double a[10][10],q[10][10],v[10][10],qt[10][10];
int i,j,n;

void cit_mat(double a[10][10]){
cout<<"\nDati dimensiunea matricei A: n=";cin>>n;
cout<<"\nDati matricea A:\n";
for (int i=1;i<=n;i++)
for (int j=1;j<=n;j++){
cout<<"a["<<i<<"]["<<j<<"]=";
cin>>a[i][j];
}
}

void afis_mat(double a[10][10],char c[1]){
cout<<"\n"<<c<<"\n";
for (int i=1;i<=n;i++){
for (int j=1;j<=n;j++){
cout.precision(4);
cout.width(4);
cout<<a[i][j]<<"   ";
}
cout<<endl;
}
}

void produs(double a[10][10],double b[10][10],double c[10][10]){
double c1[10][10];
int i,j;
for (i=1;i<=n;i++){
for (j=1;j<=n;j++){
c1[i][j]=0;
for (int k=1;k<=n;k++) c1[i][j]+=a[i][k]*b[k][j];
}
}
for (i=1;i<=n;i++)
for (j=1;j<=n;j++)
c[i][j]=c1[i][j];
}

void unitate(double a[10][10]){
for (int i=1;i<=n;i++)
for (int j=1;j<=n;j++)
if (i==j) a[i][j]=1;
else a[i][j]=0;
}

void Alg25(double a[10][10], double q[10][10]){
double h[10][10],sigma,beta,v[10],s;
int i,j,k;
unitate(q);
for (k=1;k<=n-1;k++){
s=0;
for (i=k;i<=n;i++)
s+=a[i][k]*a[i][k];
sigma=sqrt(s);
if (sigma!=0)  {
if (a[k][k]<0)
sigma=-sigma;
v[k]=a[k][k]+sigma;
for(i=1;i<=k-1;i++)
v[i]=0;
for (i=k+1;i<=n;i++)
v[i]=a[i][k];
s=0;
for(i=k;i<=n;i++)
s+=v[i]*v[i];
beta=s/2;
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
if (i!=j) h[i][j]=-v[i]*v[j]/beta;
else h[i][j]=1-v[i]*v[i]/beta;
produs(h,a,a);
produs(q,h,q);
}
}
}

void Alg40(double r[10][10],double v[10][10]){
double prod,s;
int i,j,k;
prod=1;
for(k=1;k<=n;k++)
prod*=r[k][k];
if (prod!=0){
for(j=1;j<=n-1;j++)
for(i=j+1;i<=n;i++)
v[i][j]=0;
v[1][1]=1/r[1][1];
for (k=2;k<=n;k++) {
v[k][k]=1/r[k][k];
for(i=k-1;i>=1;i--){
s=0;
for(j=i+1;j<=k;j++)
s+=r[i][j]*v[j][k];
v[i][k]=-s/r[i][i];
}
}
}
else{
cout<<"\nMatricea A nu este inversabila.";
exit(1);
}
}

void main(void){
clrscr();
cit_mat(a);
Alg25(a,q);
Alg40(a,v);//Matricea R este A
for(i=1;i<=n;i++)
for(j=1;j<=n;j++)
qt[i][j]=q[j][i];
produs(v,qt,a);
afis_mat(a,"Inversa matricei A este:\n");
getch();
}``````

658 total views, 1 views today