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 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128
|
#include<iostream>
#include<vector>
#include<cmath>
#include <iomanip>
using namespace std;
inline void pivot(vector<vector<double>> &b,vector<double>&r,int num,int n)
{
int p1,p2,g;
for (p1 = num;p1<=n-1;p1++)
{
for (p2 = p1+1;p2<=n;p2++)
{
if ( abs(b[p1][num]) < abs(b[p2][num]) )
{
swap(r[p1],r[p2]);
for(g=0;g<=n;g++)
{
swap (b[p1][g],b[p2][g]);
}
}
}
}
}//pivot
void gauss(vector<vector<double>> b,vector<double> r, vector<double> &ans)
{
int n=b.size(),p1,p2,num,k,row,g,z,j;
double f;
ans.resize(n);
n=n-1;
for(k=0;k<=n-1;k++)
{
pivot(b,r,k,n);
for(row=k;row<=n-1;row++)
{
if (b[row+1][k]==0) {break;};
f=b[k][k]/b[row+1][k];
r[row+1]=r[row+1]*f-r[k];
for(g=0;g<=n;g++)
{
b[row+1][g]=b[row+1][g]*f-b[k][g];
}//g
}//row
}//k
//back substitute
for(z=n;z>=0;z--)
{
ans[z]=r[z]/b[z][z];
for(j=n;j>=z+1;j--)
{
ans[z]=ans[z]-(b[z][j]*ans[j]/b[z][z]);
}//j
}//z
}//gauss
void matmultvect(vector<vector<double>> m1,vector<double> m2, vector<double> &ans)
{
int rows=m1.size(),r,k;
double rxc;
ans.resize(rows);
rows=rows-1;
for(r=0;r<=rows;r++)
{
rxc=0;
for(k=0;k<=rows;k++)
{
rxc=rxc+m1[r][k]*m2[k];
}//k
ans[r]=rxc;
}//r
}//matmult
void print(vector<double> v)
{
int n=v.size(),i;
string sp;
for(i=0;i<n;i++)
{
if(v[i]>=0) sp="+"; else sp="";
cout<<setprecision(15)<<sp<<v[i]<<endl;
}
}//print
int main()
{
int n;
vector<vector<double>> d
{
{1, 2, -3,9.8},
{-4, 5, 6,2},
{7, 8, 9,8},
{2, 0, -3,8}
};
vector<double> r
{
{4,5,6,7}
};
vector<double> ret;
gauss(d,r,ret);
cout<<"solution:"<<endl;
print(ret);
cout <<" "<<endl;
// multiply back to check
vector<double>ret2;
matmultvect(d,ret,ret2);
cout<<"check, should be 4,5,6,7:"<<endl;
print(ret2);
cout <<" "<<endl;
cout <<"Press return to end . . ."<<endl;
cin.get();
}
| |