Метод магнитостатики конечных элементов на треугольниках, страница 6

printf("k=%5d EPS=%10.5e a=%e norm=%e norm_z=%e\n",k,e2,a,x.norm(),z.norm());

}

e2=(*this*x-f).norm()/f_norm;

printf("EPS=%e\n",e2);

//VECTOR ff=*this*x;

//double aaa=(ff-f).norm()/f.norm();

}

//------------------------------------------------void SLAU::MSG2()

{

VECTOR r(n),z(n),h(n),h2(n),h3(n);

x.ins_el_all(0);

r=f-*this*x;

z=r;

double f_norm=f.norm();

double tmp;

double e2=r.norm()/f_norm;

double a,b;

long k;

tmp=r*r;

for(k=1;e2>EPS;k++)

{

//a

h=*this*z;//A*z

a=tmp/(h*z);

//r

r=r-a*h;

//x

x=x+a*z;

//b

b=1./tmp;

//h2=solution_x(r); //(SST)^-1*r в h2

tmp=r*r;

b*=tmp;

//z

z=r+b*z;

e2=r.norm()/f_norm;

printf("k=%5d EPS=%10.5e\n",k,e2);

}

e2=(*this*x-f).norm()/f_norm;

printf("EPS=%e\n",e2);

}

//-------------------------------------------------void SLAU::BCG()

{

VECTOR r(n),r1(n),p(n),p1(n),h(n);

x.ins_el_all(0);

r=f-*this*x;

r1=r;

p=r;

p1=r1;

double f_norm=f.norm();

double e2=r.norm()/f_norm;

double a,b;

double temp=r*r1;

for(long k=0;e2>EPS;k++)

{

h=*this*p;

a=temp/(h*p1);

x=x+a*p;

r=r-a*h;

r1=r1-a*(*this*p1);

b=1./temp;

temp=r*r1;

b*=temp;

if(fabs(b)<EPS) printf("ERROR IN BCG\n");

p=r+b*p;

p1=r1+b*p1;

e2=r.norm()/f_norm;

printf("k=%5d EPS=%10.5e\n",k,e2);

}

e2=(*this*x-f).norm()/f_norm;

printf("EPS=%e\n",e2);

}

//-------------------------------------------------void SLAU::LOS()

{

MATRIX s;

SST(&s);

VECTOR r(n),z(n),h(n),p(n),h2(n);

x.ins_el_all(0);

h=f-*this*x;

r=s.solution_x_l(h);

z=s.solution_x_u(r);

h=*this*z;

p=s.solution_x_l(h);

double f_norm=r.norm();

double tmp;

double e2=r.norm()/f_norm;

double a,b;

long k;

//h2=s.solution_x(r); //(SST)^-1*r в h2

//tmp=h2*r;

for(k=1;e2>EPS;k++)

{

//a

a=(p*r)/(p*p);

//x

x=x+a*z;

//r

r=r-a*p;

//b

h=s.solution_x_u(r);

h2=*this*h;

h=s.solution_x_l(h2);

b=-(p*h)/(p*p);

//z

h2=s.solution_x_u(r);

z=h2+b*z;

//p

p=h+b*p;

e2=r.norm()/f_norm;

/*VECTOR test=f-*this*x;

test=s.solution_x_l(test);

double eee=fabs(e2-test.norm()/f_norm);

if(eee>EPS)

{printf("ERROR eps-e\n");getch();}*/

printf("k=%5d EPS=%10.5e norm=%e\n",k,e2,x.norm());

}

e2=(*this*x-f).norm()/f_norm;

printf("EPS=%e\n",e2);

}

//-------------------------------------------------void SLAU::LOS2()

{

VECTOR r(n),z(n),h(n),p(n),h2(n);

x.ins_el_all(0);

r=f-*this*x;

z=r;

p=*this*z;

double f_norm=f.norm();

double tmp;

double e2=r.norm()/f_norm;

double a,b;

long k;

for(k=1;e2>EPS;k++)

{

//a

a=p*r/(p*p);

//x

x=x+a*z;

//r

r=r-a*p;

//b

h=*this*r;

b=-(p*h)/(p*p);

//z

z=r+b*z;

//p

p=h+b*p;

e2=r.norm()/f_norm;

//e2=(xx-x).norm()/x.norm();

printf("k=%5d EPS=%10.5e norm=%e\n",k,e2,x.norm());

}

e2=(*this*x-f).norm()/f_norm;

printf("EPS=%e\n",e2);

}

//-------------------------------------------------VECTOR MATRIX::solution_x(VECTOR f)

{

if(f.n!=n) {printf("ERROR IN SOLUTION_X\n");return NULL;}

VECTOR x=f;

long i;

//SY=F

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

{

for(long j=ig[i]-1;j<ig[i+1]-1;j++)

x.vect[i]-=gg[j]*x.vect[jg[j]-1];

x.vect[i]/=d[i];

}

//STX=Y

for(i=n-1;i>=0;i--)

{

x.vect[i]/=d[i];

for(long j=ig[i+1]-2;j>=ig[i]-1;j--)

x.vect[jg[j]-1]-=gg[j]*x.vect[i];

}

return x;

}

//--------------------------------VECTOR MATRIX::solution_x_l(VECTOR f)

{

if(f.n!=n) {printf("ERROR IN SOLUTION_X\n");return NULL;}

VECTOR x=f;

long i;

//SY=F

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

{

for(long j=ig[i]-1;j<ig[i+1]-1;j++)

x.vect[i]-=gg[j]*x.vect[jg[j]-1];

x.vect[i]/=d[i];

}

return x;

}

//--------------------------------VECTOR MATRIX::solution_x_u(VECTOR f)

{

if(f.n!=n) {printf("ERROR IN SOLUTION_X\n");return NULL;}

VECTOR x=f;

//STX=Y

for(long i=n-1;i>=0;i--)

{

x.vect[i]/=d[i];

for(long j=ig[i+1]-2;j>=ig[i]-1;j--)

x.vect[jg[j]-1]-=gg[j]*x.vect[i];

}

return x;

}

//--------------------------------void MATRIX::SST(MATRIX *my)

{

long i,j,k,ll,kk;

double sum;

my->n=n;

my->d=(double*)calloc(n,sizeof(double));

my->ig=(long*)calloc(n+1,sizeof(long));

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

my->ig[i]=ig[i];

my->jg=(long*)calloc(ig[n]-1,sizeof(long));