Witam.
Piszę program, który rozwiązuje układ równań metodą gradientów prostych.
Dla małych wielkości zadania wszystko jest ok, natomiast gdy podaję rozmiar większy niż 7 główna pętla nie chce się zatrzymać.
Co ciekawe dla rozmiaru = 5 też jest źle ale już dla6i7 jest ok.
Biblioteka, której używałem to Eigen, mam nadzieję, że ktoś tutaj się z nią spotkał.
Poniżej kod mojego programu.
Bardzo proszę o pomoc/sugestie.
#include
#include
#include
#include
using namespace std;
USING_PART_OF_NAMESPACE_EIGEN
int main(int, char *[])
{
cout << "Program rozwiazujacy uklad rownan metoda gradientow prostych" << endl << endl;
int rozmiar=8;
MatrixXd Los(rozmiar,rozmiar);
VectorXd B(rozmiar);
MatrixXd Temp(rozmiar,rozmiar);
MatrixXd A(rozmiar,rozmiar);
MatrixXd Dk(rozmiar,rozmiar);
VectorXd X(rozmiar);
VectorXd Adk(rozmiar);
RowVectorXd DkT(rozmiar);
RowVectorXd XT(rozmiar);
RowVectorXd BT(1,rozmiar);
MatrixXd Ma(1,1);
MatrixXd Mb(1,1);
VectorXd Xchol(rozmiar);
int i,j,licznik = 0,k;
double eps = 0.00001;
long double Alfa = 0;
bool Minimum = false;
double norm;
for (i=0;i
{ // Losowanie wartości macierzy
for (j=0;j
{
Los(i,j) = Los(j,i) = rand();
B(j) = rand();
}
}
cout << "Wylosowana macierz: " << endl << Los << " " <
Temp = Los.transpose(); // Uzyskiwanie Macierzy symetrycznej, dodatnio określonej
A = Los * Temp;;
cout << "Macierz wejsciowa zadania (sym - dod - okr) : " << endl << A << " " <
cout << "Wylosowany wektor B:" << endl << B << endl << endl;
cout << "Rozwiazanie Choleskym:" << endl << endl;
A.llt().solve(B,&Xchol);
cout << Xchol << endl;
cout << endl << "Licze............" << endl << endl;
while (Minimum == false)
{ //Głwna pętla
licznik++;
// cout << endl << licznik;
Dk = -1 * A * X + B;
//Dktemp = -1*A;
// Dk = prod(Dktemp,X)+b;
//cout << Dk; // Liczymy Dk - przeciw-gradient
norm = Dk.norm();
if(norm < eps) // Sprawdzanie warunku wyjścia z pętli (norma Dk < niż przyjęty epsilon
{
Minimum = true;
} else {
Adk = A*Dk; // Liczymy iloczyn A oraz Dk
DkT = Dk.transpose(); // Transponujemy Dk
BT = B.transpose(); // Transponujemy wektor B
Ma =DkT * Adk;
// Składnik a Alfy
XT = X.transpose();
Mb = XT * Adk - BT * Dk;
// Składnik b Alfy
Alfa = (-1 * Mb(0,0)) / (2 * Ma(0,0)); // Obliczenie Alfy
X = X + Alfa * Dk;
}
}
cout << "Liczba obrotow petli: " << licznik << endl << endl;
cout << "Rozwiazanie metoda gradientow prostych: " << endl;
cout << endl << X << endl << endl;
// cout << endl << Dk << endl;
// cout << endl << Ma << endl;
// cout << endl << Mb << endl;
// cout << endl << Alfa << endl;
//cout << endl << norm << endl;
cout << "Rozwiazanie Choleskym:" << endl << endl;
A.llt().solve(B,&Xchol);
cout << Xchol << endl;
system ("PAUSE");
}