1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83
|
#include<iostream>
using namespace std;
#include<vector>void
Pivotage(vector<double> & arg,int l,int m, int ligne, int column)
{
int i,j;
double pivot,coef,perm;
i=l;
while(arg[i*column+m]==0 && i<ligne)
{
i++;
}
if(i==ligne)
{
cerr <<"Pivot null et a pas lieu à une permutation" << endl;
return;
}
else
{
/*permuter les deux ligne*/
for(j=0;j<column;j++)
{
perm=arg[l*column+j];
arg[l*column+j]=arg[i*column+j];
arg[i*column+j]=perm;
}
}
pivot=arg[l*column+m];
for(i=0;i<ligne;i++)
if(i!=l)
{
coef=arg[i*column+m]/pivot;
arg[i*column+m]=0;
for(j=0;j<column;j++)
if(j!=m)
{
arg[i*column+j]=arg[i*column+j]-coef*arg[l*column+j];
}
}
coef=1/pivot;
for(j=0;j<column;j++)
arg[l*column+j]=coef*arg[l*column+j];
//cout<<"pivotage"<<endl;
}
vector<double> Inverse(vector<double>const & arg,int ligne)
{
int i,j,column=ligne*2;
vector<double> arg1(ligne*column);
//arg1.resize(ligne*column);
//double *arg1 = new double[ligne * column];
//double *inver=new double[ligne*ligne];
vector<double> inver(ligne*ligne);
//inver.resize(ligne*ligne);
//cout<<"le pb n'est pas dans l'inverse"<<endl;
for(i=0;i<ligne;i++)
for(j=0;j<ligne;j++)
arg1[i*column+j]=arg[i*ligne+j];
for(i=0;i<ligne;i++)
for(j=0;j<ligne;j++)
arg1[i*column+(j+ligne)]=0;
for(i=0;i<ligne;i++)
arg1[i*column+(i+ligne)]=1;
//cout<<"lign"<<ligne<<endl;
for(i=0;i<ligne;i++)
{
Pivotage(arg1,i,i,ligne,column);
}
for(i=0;i<ligne;i++)
for(j=0;j<ligne;j++)
inver[i*ligne+j]=arg1[i*column+(j+ligne)];
// delete[] arg1;
// cout<<"inverse"<<endl;
return(inver);
}
void main()
{
vector<double> A(n*n),B(n*n);
int n;
B=Inverse(A,n);
} |
Partager