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

795 total views, 1 views today