Pagini recente » Cod sursa (job #1898456) | Cod sursa (job #496312) | Cod sursa (job #2751213) | Cod sursa (job #2672246) | Cod sursa (job #1808087)
/*-----------------------------------------------------------
GaussJordan Algorithm (elimination with full pivoting)
Uses:
(1) solving systems of linear equations (AX=B)
(2) inverting matrices (AX=I)
(3) computing determinants of square matrices
Running time: O(n^3)
INPUT: a[][] = an nxn matrix
b[][] = an nxm matrix
OUTPUT: X = an nxm matrix (stored in b[][])
A^{-1} = an nxn matrix (stored in a[][])
returns determinant of a[][]
-----------------------------------------------------------*/
#include <bits/stdc++.h>
using namespace std;
const double EPS = 1e-10;
typedef vector<int> VI;
typedef double T;
typedef vector<T> VT;
typedef vector<VT> VVT;
T GaussJordan(VVT &a, VVT &b)
{
const int n = a.size();
const int m = b[0].size();
VI irow(n), icol(n), ipiv(n);
T det = 1;
for (int i = 0; i < n; i++)
{
int pj = -1, pk = -1;
for (int j = 0; j < n; j++)
if (!ipiv[j])
for (int k = 0; k < n; k++)
if (!ipiv[k])
if (pj == -1 || fabs(a[j][k]) > fabs(a[pj][pk])) { pj = j; pk = k; }
if (fabs(a[pj][pk]) < EPS)
{
cout<<"Impossible\n";
exit(0);
}
ipiv[pk]++;
swap(a[pj], a[pk]);
swap(b[pj], b[pk]);
if (pj != pk) det *= -1;
irow[i] = pj;
icol[i] = pk ;
T c = 1.0 / a[pk][pk];
det *= a[pk][pk];
a[pk][pk] = 1.0;
for (int p = 0; p < n; p++)
a[pk][p] *= c;
for (int p = 0; p < m; p++)
b[pk][p] *= c;
for (int p = 0; p < n; p++)
{
if (p != pk)
{
c = a[p][pk];
a[p][pk] = 0;
for (int q = 0; q < n; q++)
a[p][q] -= a[pk][q] * c;
for (int q = 0; q < m; q++)
b[p][q] -= b[pk][q] * c;
}
}
}
for (int p = n-1; p >= 0; p--)
{
if (irow[p] != icol[p])
{
for (int k = 0; k < n; k++)
swap(a[k][irow[p]], a[k][icol[p]]);
}
}
return det;
}
#define N 310
double a[N][N]={0};
double b[N]={0};
int main()
{
ios_base::sync_with_stdio(0);
cin.tie(0);
freopen("gauss.in","r",stdin);
freopen("gauss.out","w",stdout);
int n,m;
cin>>n>>m;
for(int i=0;i<n;++i) {
for(int j=0;j<m;++j)
cin>>a[i][j];
cin>>b[i];
}
int p = max(n,m);
VVT A(p), B(p);
for(int i=0;i<p;++i) {
A[i] = VT(a[i],a[i]+p);
B[i].push_back(b[i]);
}
T det = GaussJordan(A,B);
for(int i=0;i<m;++i)
cout<<fixed<<setprecision(10)<<B[i][0]<<" ";
cout<<endl;
}