[C++] Pętla nie zawsze się zatrzymuje


(Pinkibass) #1

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");

}

([alex]) #2

Sprawdź czy "norm" zawsze jest dodatnia.


(Pinkibass) #3

Dzięki [alex] za podpowiedź.

Niestety nie udało mi się jeszcze rozwiązać tego ale sprawdziłęm jak zmienia się norm

Okazało się, że dla tych rozmiarów jest wszystko ok, parametr norm przyjmuje odpowiednie wartości.

Dla rowmiarów, dla których jest źle, parametr norm cały czas ma wartość -1.#IO

Wie ktoś może co ona oznacza ?

Pozdro

Edit:

Dziękuję za wszystkie podpowiedzi.

Problem jest rozwiązany.

Okazało się, że przyczyną problemu był kompilator Dev C++

Po przejściu na MS Visual C++ wszystko jest OK.