matrix矩阵求逆 与解方程模板 留做备用 (有bug,待补充)

 //
// main.cpp
// 矩阵求逆
//
// Created by 唐 锐 on 13-6-20.
// Copyright (c) 2013年 唐 锐. All rights reserved.
// #include<iostream>
#include<algorithm>
#include<iomanip>
#include<string>
#include<sstream>
#include<cmath>
#include<vector>
using namespace std;
namespace MATRIX {
const double eps = 1e-;
template<typename T>
class matrix
{
vector<vector<T>>m;
int check()
{
if(m.empty())return ;
for(int i=;i<m.size();i++)
if(m[i].size()!=m[].size())
return ;
return ;
};
matrix wrong()
{
matrix mat;
mat.m.clear();
vector<T>vec(,-);
mat.m.push_back(vec);
return mat;
}
public:
T at(int i,int j)
{
return m[i][j];
}
vector<T>& operator [](const int i) const
{
return m[i];
}
int size()
{
return m.size();
}
istream &input(istream& in)
{
do {
m.clear();
string s;
while(getline(in,s))
{
if(s=="")break;
istringstream temp(s);
T t;
vector<T> v ;
while(temp>>t)
v.push_back(t);
m.push_back(v);
}
} while(!check());
return in;
}
ostream &output(ostream& out) const
{
for(int i=;i<m.size();i++)
{
for(int j=;j<m[i].size();j++)
{
out<<m[i][j]<<' ';
}
out<<endl;
}
return out;
}
friend inline istream& operator>>(istream& in,matrix& m)
{
return m.input(in);
}
friend inline ostream& operator<<(ostream& out,const matrix& m)
{
return m.output(out);
}
void eye(int n,int k=-)
{
if(k==-)k=n;
m.clear();
vector<T> vec(k,);
vector<vector<T>> mat(n,vec);
for(int i=;i<n&&i<k;i++) mat[i][i]=;
m=mat;
} matrix inv()
{
if(m.size()!=m[].size())return wrong();
matrix ans;
vector<vector<T>>cp=m;
ans.eye((int)m.size());
for(int i=;i<cp.size();i++)
{
if(fabs(cp[i][i])<eps)
for(int j=i+;j<cp.size();j++)
{
if(fabs(cp[j][i])>eps)
{
swap(cp[i],cp[j]);
break;
}
if(j==cp.size()-) return wrong();
}
for(int j=i+;j<cp.size();j++)
{
T k=cp[j][i]/cp[i][i];
for(int l=;l<cp[i].size();l++)
{
cp[j][l]-=k*cp[i][l];
ans.m[j][l]-=k*ans.m[i][l];
}
}
} for(int i=(int)cp.size()-;i>=;i--)
{
for(int j=i-;j>=;j--)
{
T k=cp[j][i]/cp[i][i]; for(int l=;l<cp[i].size();l++)
{
cp[j][l]-=k*cp[i][l];
ans.m[j][l]-=k*ans.m[i][l];
}
}
}
for(int i=;i<cp.size();i++)
for(int j=;j<cp[i].size();j++)
{
ans.m[i][j]/=cp[i][i];
}
return ans;
}
vector<T>solve(vector<T>v)
{
vector<T> wrong;
vector<int>turn;
for(int i=;i<v.size();i++)
turn.push_back(i);
if(m[].size()!=v.size())return wrong;
vector<T> ans(v.size(),);
vector<vector<T>>cp=m;
for(int i=;i<cp.size();i++)
{
if(fabs(cp[i][i])<eps)
for(int j=i+;j<cp.size();j++)
{
if(fabs(cp[j][i])>eps)
{
swap(cp[i],cp[j]);
swap(v[i],v[j]);
swap(turn[i],turn[j]);
break;
}
if(j==cp.size()-) return wrong;
}
for(int j=i+;j<cp.size();j++)
{
T k=cp[j][i]/cp[i][i];
for(int l=;l<cp[i].size();l++)
{
cp[j][l]-=k*cp[i][l];
}
v[j]-=k*v[i];
}
}
for(int i=(int)cp.size()-;i>=;i--)
{
for(int j=i-;j>=;j--)
{
T k=cp[j][i]/cp[i][i]; for(int l=;l<cp[i].size();l++)
{
cp[j][l]-=k*cp[i][l];
}
v[j]-=k*v[i];
}
}
for(int i=;i<v.size();i++)
{
ans[turn[i]]=v[i]/cp[i][i];
if(fabs(ans[turn[i]])<eps)
ans[turn[i]]=;
}
return ans;
} };
template <typename T>
ostream &operator<<(ostream& out,const vector<T>& v)
{
for(int i=;i<v.size();i++)
{
if(i)out<<' ';
out<<v[i];
}
return out;
} }
using namespace MATRIX;
int main()
{
matrix<double> m;
double t;
vector<double> v;
while(cin>>m)
{
for(int i=;i<m.size();i++)
{
cin>>t;
v.push_back(t);
}
vector<double>ans=m.solve(v);
cout<<ans<<endl;
}
} /* matrix 1 1 2 0 0
3 4 0 0
0 0 4 1
0 0 3 2 matrix 2 1 0 -1 2 1
3 2 -3 5 3
2 2 1 4 -2
0 4 3 3 1
1 0 8 -11 4 matrix 3 1 0 -1 2 1 0 2
1 2 -1 3 1 -1 4
2 2 1 6 2 1 6
-1 4 1 4 0 0 0
4 0 -1 21 9 9 9
2 4 4 12 5 6 11
7 -1 -4 22 7 8 18 matrix 4 4 2 -3 -1 2 1 0 0 0 0
8 6 -5 -3 6 5 0 1 0 0
4 2 -2 -1 3 2 -1 0 3 1
0 -2 1 5 -1 3 -1 1 9 4
-4 2 6 -1 6 7 -3 3 2 3
8 6 -8 5 7 17 2 6 -3 5
0 2 -1 3 -4 2 5 3 0 1
16 10 -11 -9 17 34 2 -1 2 2
4 6 2 -7 13 9 2 0 12 4
0 0 -1 8 -3 -24 -8 6 3 -1 vector 1 5 12 3 2 3 46 13 38 19 -21 */
上一篇:July 08th. 2018, Week 28th. Sunday


下一篇:彻底理解js中的闭包