Calculation of the inverse of a matrix A after performing an LR factorization
//Calculul inversei unei matrice A dupa realizarea unei factorizari LR
#include<iostream>
#include<conio.h>
#include<stdlib.h>
#include<math.h>
using namespace std;
double a[10][10],P[10][10],L[10][10],S[10][10],v[10][10],u[10][10];
int 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(7);
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 Alg20(double a[10][10],double L[10][10],double P[10][10],double S[10][10]){
double max,t,M[10][10];
int i,j,k,p,q;
unitate(P);
unitate(S);
unitate(L);
for (k=1;k<=n-1;k++){
max=0;
for (i=k;i<=n;i++)
for (j=k;j<=n;j++)
if (max<=fabs(a[i][j])){
max=fabs(a[i][j]);
p=i;
q=j;
}
if (a[p][q]!=0){
if (p!=k){
for (j=1;j<=n;j++){
t=a[k][j];
a[k][j]=a[p][j];
a[p][j]=t;
t=P[k][j];
P[k][j]=P[p][j];
P[p][j]=t;
}
}
if (q!=k){
for (i=1;i<=n;i++){
t=a[i][k];
a[i][k]=a[i][q];
a[i][q]=t;
t=S[i][k];
S[i][k]=S[i][q];
S[i][q]=t;
}
}
unitate(M);
for (i=k+1;i<=n;i++)
M[i][k]=a[i][k]/a[k][k];
if(p!=k)
for(j=1;j<=n;j++) {
t=M[k][j];
M[k][j]=M[p][j];
M[p][j]=t;
}
produs(L,M,L);
for (i=k+1;i<=n;i++){
for (j=k+1;j<=n;j++)
a[i][j]-=a[i][k]*a[k][j]/a[k][k];
a[i][k]=0;
}
}
else{
cout<<"Matricea A nu este inversabila";
getch();
exit(1);
}
}
produs(P,L,L);
}
void Alg39(double l[10][10],double u[10][10]){
double prod,s;
int i,j,k;
prod=1;
for(k=1;k<=n;k++)
prod*=l[k][k];
if (prod!=0){
for(i=1;i<=n-1;i++)
for(j=i+1;j<=n;j++)
u[i][j]=0;
for (k=1;k<=n-1;k++) {
u[k][k]=1/l[k][k];
for(i=k+1;i<=n;i++){
s=0;
for(j=k;j<=i-1;j++)
s+=l[i][j]*u[j][k];
u[i][k]=-s/l[i][i];
}
}
u[n][n]=1/l[n][n];
}
else{
cout<<"\nMatricea A nu este inversabila.";
exit(1);
}
}
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(){
// clrscr();
cit_mat(a);
Alg20(a,L,P,S);
Alg39(L,u);
Alg40(a,v);//Matricea R este A
produs(S,v,a);
produs(a,u,a);
produs(a,P,a);
afis_mat(a,"Inversa matricei A este:\n");
getch();
}
833 total views, 1 views today