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

Program received signal SIGSEGV, Segmentation fault - C++

Войти
Регистрация
Восстановить пароль
 
Рейтинг: Рейтинг темы: голосов - 21, средняя оценка - 4.62
ciao
0 / 0 / 0
Регистрация: 25.12.2008
Сообщений: 28
22.04.2009, 15:48     Program received signal SIGSEGV, Segmentation fault #1
Когда запускаю дебагер:
gdb ./preci
Вюдает следуще. Я не могу понять с чем это могет быть связано.
Код:

Код
GNU gdb 6.4.90-debian
Copyright (C) 2006 Free Software Foundation, Inc.
GDB is free software, covered by the GNU General Public License, and you are
welcome to change it and/or distribute copies of it under certain conditions.
Type "show copying" to see the conditions.
There is absolutely no warranty for GDB.  Type "show warranty" for details.
This GDB was configured as "i486-linux-gnu"...Using host libthread_db library "/lib/tls/i686/cmov/libthread_db.so.1".

(gdb) run
Starting program: /users/lsg2m/perev2763/Desktop/precipitation/preci
[Thread debugging using libthread_db enabled]
[New Thread -1211492672 (LWP 5046)]

Program received signal SIGSEGV, Segmentation fault.
[Switching to Thread -1211492672 (LWP 5046)]
0x080499fc in newt::fmin (this=0xbfb6cd50, x=@0xbfb6cd88) at newt.cpp:119
119       (func->*nrfuncv)(GetN(), x, fvec);
(gdb) bt
#0  0x080499fc in newt::fmin (this=0xbfb6cd50, x=@0xbfb6cd88) at newt.cpp:119
#1  0x0804d8ed in newt::calc (this=0xbfb6cd50, x=@0xbfb6cd88, check=0xbfb6cd84)
    at newt.cpp:130
#2  0x0805b94c in main () at preci.cpp:83
Где искать ошибку. Кто-нить с этим сталкивался?
Similar
Эксперт
41792 / 34177 / 6122
Регистрация: 12.04.2006
Сообщений: 57,940
22.04.2009, 15:48     Program received signal SIGSEGV, Segmentation fault
Посмотрите здесь:

C++ Segmentation fault :(
Signal 11 (SIGSEGV) C++
C++ Segmentation fault
C++ Linux Segmentation fault
C++ Linux segmentation fault(
C++ Segmentation fault
Segmentation fault C++
Segmentation fault C++
Segmentation fault C++
C++ Segmentation fault
C++ Segmentation fault
C++ Segmentation fault

Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
После регистрации реклама в сообщениях будет скрыта и будут доступны все возможности форума.
Mecid
 Аватар для Mecid
678 / 227 / 4
Регистрация: 15.10.2007
Сообщений: 1,247
22.04.2009, 16:02     Program received signal SIGSEGV, Segmentation fault #2
Может дашь код поглядеть
ciao
0 / 0 / 0
Регистрация: 25.12.2008
Сообщений: 28
22.04.2009, 16:12  [ТС]     Program received signal SIGSEGV, Segmentation fault #3
Код
#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);

	}

}
Yandex
Объявления
22.04.2009, 16:12     Program received signal SIGSEGV, Segmentation fault
Ответ Создать тему
Опции темы

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