Dzień dobry,
obecnie zajmują się przeciążaniem operatora "/" w celu łatwego napisania mnożenia macierzy A przez odwrotność macierzy B (nieformalnie: dzielenie macierzy A przez macierz B). Przetestowałem działanie funkcji na odwracalnej macierzy, pisząc Ab/Ab, więc powinienem uzyskać macierz jednostkową. Część obliczeń się zgadza, ale w pewnym etapie, np. C[1][0], element przyjmuje pewną dziwną wartość i psuje wynik. Myślałem, że to przez przybliżenia, ale te same działania wpisane w strumień wyjścia dają wynik prawidłowy, czyli w tym przypadku 0.
Uprzejmie proszę o pomoc.
#include <iostream>
#include <vector>
#include <complex>
using namespace std;
void createIdentifyMatrix(vector<vector<complex<long double>>> &A, const unsigned n)
{
if(n==0)
{
cerr << "Wymiar takiej macierzy nie moze byc rowny 0." << endl << endl;
return;
}
const vector<complex<long double>> b(n);
for(unsigned i=0; i<n; i++)
{
A.push_back(b);
A[i][i]=1;
}
}
void printMatrix(const vector<vector<complex<long double>>> &a)
{
unsigned i, j;
for(i=0; i<a.size(); i++)
{
for(j=0; j<a[i].size(); j++)
{
cout << a[i][j] << " ";
}
cout << endl;
}
}
void matrixInversion(vector<vector<complex<long double>>> &A)
{
const unsigned n=A.size();
vector<vector<complex<long double>>> I;
complex<long double> num;
unsigned i, j, k;
createIdentifyMatrix(I, n);
for(j=0; j<n; j++)
{
k=j;
for(i=j+1; i<A.size(); i++)
{
if(fabs(A[i][j])>fabs(A[k][j]))
{
k=i;
}
}
if(A[k][j]==complex<long double>(0, 0))
{
cout << "Uklad osobliwy." << endl;
throw;
}
if(A[j][j]!=A[k][j])
{
swap(A[j], A[k]);
swap(I[j], I[k]);
}
num=A[j][j];
A[j][j]=1;
for(i=0; i<n; i++)
{
I[j][i]/=num;
if(i==j) continue;
A[j][i]/=num;
}
for(i=0; i<n; i++)
{
if(i==j) continue;
num=A[i][j]/A[j][j];
A[i][j]=0;
for(k=0; k<=j; k++)
{
I[i][k]-=I[j][k]*num;
}
for(k=j+1; k<n; k++)
{
A[i][k]-=A[j][k]*num;
I[i][k]-=I[j][k]*num;
}
}
}
A=I;
}
vector<vector<complex<long double>>> operator /(const vector<vector<complex<long double>>> &A, const vector<vector<complex<long double>>> &B)
{
const unsigned r=A.size(), c=B.size();
vector<vector<complex<long double>>> C(r, vector<complex<long double>>(c)), B1=B;
matrixInversion(B1);
unsigned i, j, k;
for(i=0; i<r; i++)
{
for(j=0; j<c; j++)
{
cout << "C[" << i << "][" << j << "]" << endl;
for(k=0; k<c; k++)
{
cout << C[i][j] << "+";
C[i][j]+=A[k][i]*B1[j][k];
cout << A[k][i] << "*" << B1[j][k] << "=" << C[i][j] << endl;
}
cout << endl;
}
}
return C;
}
complex<long double> c(const long double Re, const long double Im=0) {return complex<long double>(Re, Im);}
int main()
{
const vector<vector<vector<complex<long double>>>>
G{
{
{0,1,2,3,5},
{1,0,1,2,4},
{1,2,3,4,6},
{0,1,3,2,4}
},
{
{2,4,2,0,4},
{1,0,-1,1,2},
{0,1,3,-1,0},
{2,1,2,1,6}
},
{
{1,2,-1,2},
{2,1,0,3},
{-1,1,2,2}
},
{
{1,2,3,1},
{4,5,6,3},
{7,8,9,2}
}
},
Ch{
{
{1,1,-2,-2},
{1,2,-1,-4},
{-2,-1,6,3},
{-2,-4,3,11}
},
{
{4,2,-2},
{2,10,-7},
{-2,-7,6}
},
{
{9,-9,-6,9},
{-9,13,10,-11},
{-6,10,17,-5},
{9,-11,-5,15}
},
{
{c(4, 0), c(0, -2), c(2, -2)},
{c(0, 2), c(17, 0), c(9, 5)},
{c(2, 2), c(9, -5), c(8, 0)}
}
};
vector<vector<complex<long double>>> Ab=Ch[1];
printMatrix(Ab/Ab);
return 0;
}
