Cod sursa(job #1808080)

Utilizator afkidStancioiu Nicu Razvan afkid Data 17 noiembrie 2016 12:12:09
Problema Algoritmul lui Gauss Scor 10
Compilator cpp Status done
Runda Arhiva educationala Marime 3.22 kb
/*-----------------------------------------------------------
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; }               
        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);
    if(det==0) {
        cout<<"Imposibil\n";
        return 0;
    }
    for(int i=0;i<B.size();++i)
        cout<<fixed<<setprecision(10)<<B[i][0]<<" ";
    cout<<"\n";
    return 0;
}