Форум программистов, компьютерный форум CyberForum.ru

Program received signal SIGSEGV, Segmentation fault - C++

Войти
Регистрация
Восстановить пароль
Другие темы раздела
C++ Тестирующие сайты.. http://www.cyberforum.ru/cpp/thread31034.html
Кто решает задачи на этих тестирующих сайтах, отзовитесь и помогите решить некоторые задачи:, www.********, www.acmu.ru
C++ Комбинаторика... Обмен валюты (Время: 1 сек. Память: 16 Мб Сложность: 70%) Петя работает в обменном пункте во Флатландии. Недавно Петя получил от начальства набор цифр для отображения обменного курса. К сожалению, набор содержит всего по две копии каждой цифры. Теперь Петя хочет узнать, сколько различных обменных курсов он сможет отобразить. Петя обменивает флатландские доллары на крайландские тугрики. Петя... http://www.cyberforum.ru/cpp/thread31033.html
Использование файлов и строк C++
Уважаемые программисты! Нужна ваша помощь, заканчиваю написание диплома и возникла следующая проблемка: Есть программа, которая все результаты записывает в файл txt. Этот файл имеет следующий вид: time 0.000 0.000 0.250 0.000 0.000 0.000 0.000 ... time 0.050 0.000 0.250 0.000 0.000 0.000 0.000 ... time 0.100 0.000 0.250 0.000 0.000 0.000 0.000 ... ...
C++ Составить функцию конкатенации (слияния) двух строк
составить функцию конкатенации (слияния) двух строк.в основной програме использовать эту функцию для слияния четырех строк.
C++ Структура самоадресуемая http://www.cyberforum.ru/cpp/thread30885.html
Приветствую всех! Есть задача:чтение двоичного файла с зараннее известной структурой. Особенность:структура самоадресуемая(динамическое выделение памяти). Кто-нибудь сталкивался с такой проблемой?Инет по этому поводу молчит,Архангельский только выдает примеры на создание двоичного файла,а по поводу чтения неособо. //самоадресуемая структура struct TEtl { AnsiString Meropr; AnsiString...
C++ Какой библиотеки не хватает для работы функции sound() #include<conio.h> // какой библиотеки нехватает, чтобы интерпретатор не указывал ошибку на sound void interrupt (*SvInt09)(void); /* старый обработчик */ void interrupt NewInt09(void) { sound(100); delay(5); nosound(); SvInt09(); /* вызвать старый обработчик BIOS */ } ... void main(void) подробнее

Показать сообщение отдельно
ciao
0 / 0 / 0
Регистрация: 25.12.2008
Сообщений: 28
22.04.2009, 16:12  [ТС]     Program received signal SIGSEGV, Segmentation fault
Код
#include "erreurMatrice.h"
#ifndef newtH

#define newtH

#define ALFnewt 1.0e-4

#define epsilonnewt 1.0e-10



#define TOLXnewt epsilonnewt

#define FMAXnewt std::max

#define SQRnewt(a) ((a)*(a))

#define STPMXnewt 1.0


class newt

{
private:
  double fmin(std::vector <double> &x);
  void lnsrch(int n, std::vector <double> &xold, double fold, std::vector <double> &g, std::vector <double> &p, std::vector <double> &x,double *f, double stpmax, int *check); 
  void echange(double &a,double &b) const;
public:
  interface *func;
  interface *jacobian;
  sortie *print;
  void (interface::*nrfuncv)(int nn, std::vector <double> &v, std::vector <double> &f);
  void (interface::*jacobfunc)(std::vector <double> &, std::vector <double> &, std::vector < std::vector<double> > &);

  void (sortie::*on_subproduct)();//Function of a conclusion of intermediate results
public:
  int nn;
  std::vector <double> fvec;
  newt() //Konstruktor
  {
    on_subproduct = NULL; //By default intermediate results are not deduced
  }

  unsigned int GetN() const {return nn;};
  double norme2(std::vector <double> &V);
  void SwapRows(int i1, int i2,std::vector < std::vector<double> > &A);
  void IdentityMatrix(int n, std::vector < std::vector<double> > &A);
  void LUP(std::vector < std::vector<double> > &C, std::vector < std::vector<double> > &P);
  void Inverse(std::vector < std::vector<double> > &X);
  void calc(std::vector <double> &x, int *check); //Newton
};

//---------------------------------------------------------------------------

#endif
Код
#define BZ_THREADSAFE 
#include <blitz/array.h>


using namespace blitz;

class Gibbs;            //thermodynamyque
class interface;

class distribution;
template <class Type> class distribution_moyenne;
template <class Type> class distribution_euler;

class phase;
class matrice;
class sphere;

class siteprecipitation;
class site_homogene_constant;


class siteprecipitation
{
public:
  siteprecipitation();
  siteprecipitation(phase *);

  phase *precipiteparent;
  double densite_initiale[10];
  double densite_disponible[10];

  virtual void evolution(double)=0;
};

class site_homogene_constant : public siteprecipitation
{
public:
  site_homogene_constant();
  site_homogene_constant(phase *);

  virtual void evolution(double);
};

/*
class grain : public siteprecipitation
{
  //constructeurs
  double taillegrain[10];
  jointdegrain *jdg[10];
};


class jointdegrain : public grain
{//constructeurs
  double surface;
  double energie;
  double desorientation;
};

*/

class distribution
{
public:
  distribution();
  distribution(siteprecipitation *);
  
  phase *phasemere;
  siteprecipitation *siteparent;

  double fluxgermination;

  virtual void evolution(double)=0;
  virtual void germination_d(double)=0;
};


template <class Type>
class distribution_moyenne : public distribution
{
public:
  distribution_moyenne();
  distribution_moyenne(siteprecipitation *);  

  Type *classe0;
  
  virtual void evolution(double);
  virtual void germination_d(double);
};


template <class Type>
class distribution_euler : public distribution
{
public:
  distribution_euler();
  distribution_euler(siteprecipitation *);
  double deltaR;
  double Rcritique;
  
  Type *classe[10];
  
  virtual void evolution(double);
  virtual void germination_d(double);
};


class interface
{
public:
  interface();
  interface(phase *,phase *); //Constructeur
  phase *ph1;
  phase *ph2;
  double xeq1[10];
  double xeq2[10];
  double param_cinetique;
  void equilibre_Fe_Cr_C(int, std::vector <double> &Y, std::vector <double> &F);
  void Jacobian(std::vector <double> &Y, std::vector <double> &F, std::vector < std::vector<double> > &Jacobian2D);
};


class phase
{
public:
  phase();
  phase(distribution *);

  siteprecipitation *site[10];
  distribution *distri[10];
  phase *produit[10];

  phase *voisin[10];

  phase *phasemere;
  distribution *distributionmere;

  interface *frontiere[10];
  Gibbs *Thermo;

  double xintm[10];	//Concentration d'equilibre a l'interface dans la matrice
  double xintp[10];     //Concentration d'equilibre a l'interface dans le precipite
  double xini[10];	//Concentration initiale
  double xmoy[10];	//Concentration moyenne
  double omega[10];	//Sursaturation
  double D0[10];	//Coefficient de diffusion
  double E0[10];
  double D[10];
  double A[10];         //Coefficient cinetique
  double (*potentiel[10]) (double,double); //potentiels chimiques

  double Taille[10];
  double dTaille[10];

  double fraction_apparente;	//densite apparente * volume individuel des precipites
  double fraction;	//Fraction volumique reelle
  double volmolaire;
  double parametre_maille;
  double diffusion;
  double module_Young;
  double coefficient_Poisson;
  
  void ini(double);
  void calc_diffusion(double);
  
  virtual void calc_composition()=0;
  virtual void calc_fraction()=0;
  virtual void evolution(double)=0;
  virtual void calc_germination(double)=0;
  virtual void calc_croissance(double)=0;
};


class matrice : public phase
{
public:
  //Constructeur
  matrice();
  //Fonctions
  virtual void calc_composition();
  virtual void calc_fraction();
  virtual void evolution(double);
  virtual void calc_germination(double);
  virtual void calc_croissance(double);
};


class sphere : public phase
{
public:
  //Constructeurs
  sphere();
  sphere(distribution *);

  //Donnee d'entree
  double eninterface;
  double K0;//Constante d'equilibre
  double Q;
  double Keq;
  //Parametres microstructuraux
  //  double rayon;
  double densite_apparente;
  double densite;
  //Parametre intermediaire
  double effet_Gibbs_Thomson;		
  //Parametre d'evolution
  // double drayon;//vitesse de croissance en m.s-1
  //Fonctions
  virtual void calc_composition();
  virtual void calc_fraction();		
  virtual void evolution(double);
  virtual void calc_germination(double);		
  virtual void calc_croissance(double);
  void copie(sphere *);
};


class sphere_voisin : public sphere
{
public:
  //Constructeurs
  sphere_voisin();
  sphere_voisin(phase *);

  double drayon2;
  double rayon2;

  virtual void evolution(double);
  virtual void calc_croissance(double);
};

class sphere_germ_homogene : public sphere
{
public:
  //Constructeurs
  sphere_germ_homogene();
  sphere_germ_homogene(distribution *);

  //Donnees d'entree
  double eninterfacegerm;
  double vatomique;
  //Parametre de microstructure
  double nsites;
  //Parametres intermediaires pour le calcul du flux de germination
  double forcemotrice;//par unite de volume
  double zeldovich;
  double transport;
  double barriere;
  double rayoncritique;
  //Parametre d'evolution
  double fluxgermination;
  //Fonctions		
  virtual void evolution(double);
  virtual void calc_germination(double);
};

class sphere_germ_homogene_elasticite : public sphere_germ_homogene
{
public:
  //Constructeurs
  sphere_germ_homogene_elasticite();
  sphere_germ_homogene_elasticite(distribution *);

  //Donnee d'entree
  double misfit;
  //Parametre intermediaire
  double effet_energie_elastique;
  //Fonctions
  virtual void calc_germination(double);
};

//=====================================Class Gibbse=========================================================
class Gibbs
{  
FILE *fp;
double T1, T2;
private:
double polynom (double T1,Array <double,1> &a);
public:
  Gibbs();
  Gibbs(const char *);
  const char *datafile;

  std::vector <double> y1;
  std::vector <double> y2;
  std::vector <double> m; 
  double Gm;
  std::vector <double> dGm1, dGm2;
  std::vector <double> pot_chemical1, pot_chemical2;
  double Tc, Bmag; 
//===================================================================

 int calc_spaces(char *);
 void read(double);           

//===================================================================
void rewriteGHSERp(double,Array <double,1> &, Array <double,4> &, Array <double,3> &); 
void rewriteL(double, Array <double,1> &, Array <double,6> &, Array <double,5> &);
void rewriteLtetra(double, Array <double,1> &, Array <double,7> &, Array <double,6> &);

void rewriteGHSERfinal(double, Array <double,3> &, Array <double,2> &);
void rewriteLfinal(double, Array <double,5> &, Array <double,4> &);
void rewriteL_tetra_final(double, Array <double,6> &,Array <double,5> &);

void scanmatrixL(Array <double,4> &, Array <double,1> &, Array <int,1> &, Array <int,1> &, Array <int,1> &, Array <int,1> &);
void scanmatrixLtetra(Array <double,5> &,Array <double,1> &,Array <int,1> &,Array <int,1> &,Array <int,1> &,Array <int,1> &,Array <int,1> &);

//===================================================================
void functionGibbs(double,Array <double,2> &,Array <double,1> &,Array <int,1> &,Array <int,1> &,Array <int,1> &,Array <int,1> &,Array <double,1> &,Array <int,1> &,Array <int,1> &,Array <int,1> &,Array <int,1> &,Array <double,1> &,Array <int,1> &,Array <int,1> &,Array <int,1> &,Array <int,1> &, Array <int,1> &);

void function_derivation_Gibbs(double T,Array <double,2> &Gfinal,Array <double,1> &L1,Array <int,1> &iL1,Array <int,1> &jL1,Array <int,1> &kL1,Array <int,1> &lL1,Array <double,1> &L2,Array <int,1> &iL2,Array <int,1> &jL2,Array <int,1> &kL2,Array <int,1> &lL2,Array <double,1> &Ltetra,Array <int,1> &iLtetra,Array <int,1> &jLtetra,Array <int,1> &kLtetra,Array <int,1> &lLtetra, Array <int,1> &mLtetra,Array <double,2> &TCfinal,Array <double,1> &TC1,Array <int,1> &iTC1,Array <int,1> &jTC1,Array <int,1> &kTC1,Array <int,1> &lTC1,Array <double,1> &TC2,Array <int,1> &iTC2,Array <int,1> &jTC2,Array <int,1> &kTC2,Array <int,1> &lTC2,Array <double,1> &TCtetra,Array <int,1> &iTCtetra,Array <int,1> &jTCtetra,Array <int,1> &kTCtetra,Array <int,1> &lTCtetra, Array <int,1> &mTCtetra,Array <double,2> &BMAGfinal,Array <double,1> &BMAG1,Array <int,1> &iBMAG1,Array <int,1> &jBMAG1,Array <int,1> &kBAMG1,Array <int,1> &lBMAG1,Array <double,1> &BMAG2,Array <int,1> &iBMAG2,Array <int,1> &jBMAG2,Array <int,1> &kBMAG2,Array <int,1> &lBMAG2,Array <double,1> &BMAGtetra,Array <int,1> &iBMAGtetra,Array <int,1> &jBMAGtetra,Array <int,1> &kBMAGtetra,Array <int,1> &lBMAGtetra, Array <int,1> &mBMAGtetra);

void chemical_potential();

double functionTC_BMAG(double,Array <double,2> &,Array <double,1> &,Array <int,1> &,Array <int,1> &,Array <int,1> &,Array <int,1> &,Array <double,1> &,Array <int,1> &,Array <int,1> &,Array <int,1> &,Array <int,1> &,Array <double,1> &,Array <int,1> &,Array <int,1> &,Array <int,1> &,Array <int,1> &, Array <int,1> &);

};
Добавлено через 1 минуту 13 секунд
Код
#include "donnees.cpp"
#include "phase.cpp"
#include <ncurses.h>
#include "newt.cpp"
#include "erreurMatrice.cpp"


char fichiersortie[20];
segment seg[20];
int nbseg;
chargement ch;
sortie so;

matrice mat;

site_homogene_constant site_cem(&mat);
site_homogene_constant site_eps(&mat);

distribution_moyenne<sphere> cem(&site_cem);
distribution_moyenne<sphere> eps(&site_eps);

double temps,dt;
double T,dT;

int main()
{
mat.Thermo = new Gibbs("ferrite.txt");//Fe-Cr-C
cem.classe0->Thermo = new Gibbs("cementite.txt");//Fe-Cr-C
newt nt;
phase *matrice;
interface mv;
interface pf;


  ch.ouverture("param.txt");
  ch.lecture();
  ch.fermeture();
  so.ouverture(fichiersortie);//so.ouverture("sortie.txt");

    for(int iterseg=1;iterseg<=nbseg;iterseg++)
      {
        temps=seg[iterseg].date;
        dt=(seg[iterseg+1].date-temps)/seg[iterseg+1].pas;
      
        T=seg[iterseg].temperature;
        dT= dt * (seg[iterseg+1].temperature-seg[iterseg].temperature)
	/ (seg[iterseg+1].date-seg[iterseg].date);
      
        while(temps<seg[iterseg+1].date)
	 {
	  temps+=dt;
	  T+=dT;
//===========Newton-Raphson======
        nt.nn = 9;
        std::vector <double> Y(nt.nn);
        nt.fvec.resize(nt.nn);
        int check = 0;
        int i=1;
        while(matrice->frontiere[i]!=NULL)
         {
          nt.func = &mv;
          nt.nrfuncv = &interface::equilibre_Fe_Cr_C;
          nt.jacobian = &pf;
          nt.jacobfunc = &interface::Jacobian;
          i++;
         }
//==Thermo====
        mat.Thermo->read(T);
        cem.classe0->Thermo->read(T);
        Y[0] = mat.Thermo->y1[0];
        Y[1] = mat.Thermo->y1[1];
        Y[2] = mat.Thermo->y2[0];
        Y[3] = mat.Thermo->y2[1];
        Y[4] = cem.classe0->Thermo->y1[0];
        Y[5] = cem.classe0->Thermo->y1[1];
        Y[6] = cem.classe0->Thermo->y2[0];
        Y[7] = cem.classe0->Thermo->y2[1];
        Y[8] = nf;
//============
//        so.ecriture();
        nt.print = &so;
        nt.on_subproduct = &sortie::ecriture; 
        nt.calc(Y, &check);   
	 } 
      }
  so.fermeture();

  return(0);
}


void chargement::lecture()
{
  sphere *tempcem,*tempeps;
  tempcem=cem.classe0;
  tempeps=eps.classe0;
  
  char buf[80];
  lignesuivante();lignesuivante();lignesuivante();
  
  fscanf(fe,"%lf\n",&mat.D0[2]);
  lignesuivante();
  fscanf(fe,"%lf\n",&mat.E0[2]);
  lignesuivante();
  fscanf(fe,"%lf\n",&mat.xini[2]);
  mat.xmoy[2]=mat.xini[2];
  lignesuivante();
  fscanf(fe,"%lf\n",&mat.volmolaire);
  lignesuivante();	
  fscanf(fe,"%lf\n",&site_cem.densite_initiale[1]);
  cem.classe0->densite_apparente=site_cem.densite_initiale[1];
  site_cem.densite_disponible[1]=1.;
  lignesuivante();
  fscanf(fe,"%lf\n",&tempcem->Taille[1]);
  lignesuivante();
  fscanf(fe,"%lf\n",&tempcem->xmoy[2]);
  lignesuivante();
  fscanf(fe,"%lf %lf\n",&tempcem->K0,&tempcem->Q);
  lignesuivante();
  fscanf(fe,"%lf\n",&tempcem->eninterface);
  lignesuivante();
  fscanf(fe,"%lf\n",&tempcem->volmolaire);
  lignesuivante();
  fscanf(fe,"%i\n",&nbseg);
  lignesuivante();
  for(int i=1;i<=1+nbseg;i++)
    {
      fscanf(fe,"%lf %lf %lf\n",&seg[i].date,&seg[i].temperature,&seg[i].pas);
    }
  lignesuivante();lignesuivante();
  fscanf(fe,"%s\n",fichiersortie);

	
  //	tempeps->densite_apparente=tempcem->densite_apparente*2.;
	site_eps.densite_initiale[1]=site_cem.densite_initiale[1]*2.;
	eps.classe0->densite_apparente=site_eps.densite_initiale[1];
	site_eps.densite_disponible[1]=1.;
	tempeps->Taille[1]=1.e-8;
	tempeps->xmoy[2]=tempcem->xmoy[2]*2.;
	tempeps->eninterface=tempcem->eninterface*2.;
	tempeps->volmolaire=tempcem->volmolaire*2.;
	tempeps->potentiel[1]=tempcem->potentiel[1];
	tempeps->potentiel[2]=tempcem->potentiel[2];
}


void sortie::ecriture()
{
  fprintf(fs,"%i %e %e %e %e %e %e\n",
	  lignenb,//----------------------------------------1
	  temps,//------------------------------------------2
	  T,//----------------------------------------------3
          cem.classe0->Thermo->pot_chemical1[0],//----------4
          cem.classe0->Thermo->pot_chemical1[1],//----------5
          cem.classe0->Thermo->pot_chemical2[0],//----------6
          cem.classe0->Thermo->pot_chemical2[1],//----------7
          mat.Thermo->pot_chemical1[0],//-------------------8
          mat.Thermo->pot_chemical1[1],//-------------------9
          mat.Thermo->pot_chemical2[0],//-------------------10
          mat.Thermo->pot_chemical2[1]//-------------------11
);   
lignenb+=1;
}
Добавлено через 2 минуты 0 секунд
Выслала только структуру проги. А так она очень большая.
С чем проблемма связана, не могу понять.

Добавлено через 2 минуты 30 секунд
Работаю с классом interface i Gibbs.
Уравнения, которие долген решать метод Ньютон-Рафсон записываю в класс interface.

Добавлено через 1 минуту 7 секунд
Код
#include "newt.h"

void newt::echange(double &a,double &b) const

{
   double tmp = a;

   a = b;

   b = tmp;
}

//========================Norm of a vector=========================================================

double newt::norme2(std::vector <double> &V)

{

	double norme = 0.0;

	for(unsigned int i=0; i<nn; i++)

		norme += V[i]*V[i];



	return sqrt(norme);

}

//==============Transpose of a matrice (Lines are interchanged the position i1 and i2)======================================
void newt::SwapRows(int i1, int i2,std::vector < std::vector<double> > &A)

{

  assert((i1 < nn) && (i2 < nn));

  for (int j = 0; j < nn; ++j)

    echange(A[i1][j], A[i2][j]);

}
//==============Return of an identical matrix===========================================================
void newt::IdentityMatrix(int n, std::vector < std::vector<double> > &A)
{
  std::vector < std::vector<double> > res(n);
  for(int i=0 ; i < res.size() ; i++)
  {
   res[i].resize(n);
  }
  for(int i=0 ; i < res.size() ; i++)
  res[i][i] = 1;
  A = res;
}

//==============LUP-decomposition (C= L + U - E)=======================================================
void newt::LUP(std::vector < std::vector<double> > &C, std::vector < std::vector<double> > &P)
{
  IdentityMatrix(nn, P);

  for( int i = 0; i < nn; i++ ) {
  double pivotValue = 0;
  int pivot = -1;
  for( int row = i; row < nn; row++ ) {
    if( fabs(C[row][i]) > pivotValue ) {
       pivotValue = fabs(C[row][i]);
       pivot = row;
         }
     }
    if( pivotValue == 0 ) {

    cout<<"matrix is singular"<<endl;
    return;
        }

    SwapRows(pivot, i, P);
    SwapRows(pivot, i, C);
    for( int j = i+1; j < nn; j++ ) {
    C[j][i] /= C[i][i];
      for( int k = i+1; k < nn; k++ )
        C[j][k] -= C[j][i] * C[i][k];
        }
    }

}

//==============Inverse matrix a method LUP-decomposition=================================================
void newt::Inverse(std::vector < std::vector<double> > &X)
{
  X.resize(nn);
  for(int i = 0 ; i < X.size() ; i++)
   {
    X[i].resize(nn);
   } 
  std::vector < std::vector<double> > P(nn);
  std::vector < std::vector<double> > C(nn);
  for(int i = 0; i < P.size(); i++)
   {
    P[i].resize(nn);
   }
  for(int i = 0; i < C.size(); i++)
   {
    C[i].resize(nn);
   }
  LUP(C, P);

   for(int k = nn-1; k >= 0; k--) {

        X[k][k] = 1;

        for( int j = nn-1; j > k; j--) X[k][k] -= C[k][j] * X[j][k];

        X[k][k] /= C[k][k];

        for( int i = k-1; i >= 0; i-- ) {

            for( int j = nn-1; j > i; j-- ) {

                X[i][k] -= C[i][j] * X[j][k];

                X[k][i] -= C[j][i] * X[k][j];

            }

            X[i][k] /= C[i][i];

        }

    }

if (int m = X.size())
{

 for(int i=0 ; i < X.size() ; i++)
  for(int j=0 ; j < X[i].size() ; j++)
   {
   for (int k = 0; k < m; k++)  
     { X[i][j] += X[i][k] * P[k][j];}
   }
} 
 

}


double newt::fmin(std::vector <double> &x)
{

  int i;

  double sum;

  (func->*nrfuncv)(GetN(), x, fvec);

  for (sum = 0.0, i = 0; i < nn; ++i)

    sum += fvec[i]*fvec[i];

  return 0.5*sum;

}


void newt::calc(std::vector <double> &x, int *check)

{
  FILE *fp = fopen("sortie.txt","w");
  std::vector <double> g(nn), xold(nn);
  fvec.resize(nn);

  double fold, F = fmin(x);

  double sum = 0;

  for (int i = 0; i < nn; ++i)

  sum += SQRnewt(x[i]);

  double stpmax = STPMXnewt*FMAXnewt(sqrt(sum), (double)nn);
  std::vector <double> f(nn);
  try

  {

    while (true)

    {
      (func->*nrfuncv)(nn, x, f);
      std::vector < std::vector<double> > J(nn);
      for(int i = 0 ; i < J.size() ; i++)
      {
       J[i].resize(nn);
      } 
      (jacobian->*jacobfunc)(x,f,J);


      for (int i = 0; i < J.size(); ++i)
      {
        sum = 0;
        for (int j=0; j < J[i].size(); ++j)
          sum +=J[j][i]*fvec[j];
        g[i] = sum;
      }
      Inverse(J); //Inverse matrice Jacobian
      std::vector <double> dx(nn);

      int indice = 0;
      for (int k = 0; k < nn; ++k)
       {
        dx[indice]=0.0;
        for (int l = 0; l < nn; ++l)
           dx[l] += (-1.) * J[k][l] * f[l];
           indice++; 
       }
      fold = F;
      for (int i = 0; i < nn; ++i)
      xold[i] = x[i];
      lnsrch(nn, xold, fold, g, dx, x, &F, stpmax, check);

/************************************************************************************************************/
    if (on_subproduct) //If it is necessary to deduce results 
      (print->*on_subproduct)(); //it is deduced


     double test = 0;
     for (int i = 0; i < nn; ++i)
       if (fabs(fvec[i]) > test)
         test = fabs(fvec[i]);
     if (test < epsilonnewt)
       break;

     double norm = norme2(dx); //We consider norm

     if (norm < epsilonnewt) //If < eps 
       break; //We leave
    }

  }

  catch (ErreurMatrice &a)

  {

    printf("Errore!!!"); //If there was an error

    *check = 1;

    return;

  }
  fclose(fp);

  *check = 0;

}

void newt::lnsrch(int n, std::vector <double> &xold, double fold, std::vector <double> &g, std::vector <double> &p, std::vector <double> &x,double *f, double stpmax, int *check)

{
  	int i;

	double a,alam,alam2,alamin,b,disc,f2,fold2,rhs1,rhs2,slope,sum,temp,

		test,tmplam;



	*check=0;

	for (sum=0.0, i=0; i<n; i++) sum += p[i]*p[i];

	sum=sqrt(sum);

	if (sum > stpmax)

		for (i=0; i<n; i++)  p[i] *= stpmax/sum;
	for (slope=0.0, i=0; i<n; i++)

		slope += g[i] * p[i];

	test=0.0;

	for (i=0; i<n; i++) {

		temp = fabs(p[i])/FMAXnewt(fabs(xold[i]),1.0);

		if (temp > test) test = temp;

	}

	alamin = TOLXnewt/test;

	alam = 1.0;

	for (;;) {

		for (i=0; i<n; i++) x[i] = xold[i] + alam * p[i];

                *f = fmin(x);

		if (alam < alamin) {

//			for (i=0;i<n;i++) x[i]=xold[i];

			*check=1;

			return;

		} else if (*f <= fold+ALFnewt*alam*slope) return;

		else {

			if (alam == 1.0)

				tmplam = -slope/(2.0*(*f - fold - slope));

			else {

				rhs1 = *f - fold - alam * slope;

				rhs2=f2 - fold2 - alam2 * slope;

				a=(rhs1/(alam * alam) - rhs2/(alam2 * alam2))/(alam - alam2);

				b=(-alam2 * rhs1/(alam * alam) + alam * rhs2/(alam2 * alam2))/(alam - alam2);

				if (a == 0.0) tmplam = -slope/(2.0 * b);

				else {

					disc=b * b - 3.0 * a * slope;

					if (disc<0.0)

                                          printf("Roundoff problem in lnsrch.\n");

					else

                                          tmplam=(-b + sqrt(disc))/(3.0 * a);

				}

				if (tmplam > 0.5 * alam)

					tmplam=0.5 * alam;

			}

		}

		alam2 = alam;

		f2 = *f;

		fold2 = fold;

		alam = FMAXnewt(tmplam,0.1*alam);

	}

}
 
Текущее время: 11:14. Часовой пояс GMT +3.
КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2017, vBulletin Solutions, Inc.
Рейтинг@Mail.ru