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();
}

791 total views, 1 views today