#define MAX 2
#define ET 1e-4
#include<stdio.h>
#include<math.h>
 
main()
{
FILE *fp;
int i,j,k,count=0;
double a[MAX][MAX],ani[MAX][MAX],b[MAX],x[MAX],g[MAX],gtemp1[MAX]={0},gtemp2[MAX],ftemp[MAX],gtemp[MAX];
double p,q,s,t,f,w=1,c,lmd;
fp=fopen("d:/mnewton.txt","w+");

clrscr();
printf("Please input Hesse juzhen a:\n");/*輸入Hesse矩陣*/
for(i=0;i<MAX;i++)
  for(j=0;j<MAX;j++)
  scanf("%lf",&a[i][j]);

/*為減少程序代碼，這里略去Hesse矩陣求逆過程，由人工算出*/
printf("Please input hesse juzhen de ni juzhen:\n");
for(i=0;i<MAX;i++)
  for(j=0;j<MAX;j++)
  scanf("%lf",&ani[i][j]);

printf("Please input  juzhen b:\n");/*輸入一次項系數矩陣b*/
for(i=0;i<MAX;i++)
  scanf("%lf",&b[i]);

printf("Please input  changliang c:\n");/*輸入常數c*/
scanf("%lf",&c);

printf("please input the first point x:\n");/*輸入初始節點x*/
for(i=0;i<MAX;i++)
  scanf("%lf",&x[i]);


while(w>ET)                /*共軛梯度法的控制循環*/
{

 for(i=0;i<MAX;i++)  
 {
 ftemp[i]=0;
 gtemp1[i]=0;
 gtemp2[i]=0;
 gtemp[i]=0;
 }
 f=0;
 p=0;
 q=0;
 s=0;

 for(i=0;i<MAX;i++)   
   for(j=0;j<MAX;j++)
   ftemp[i]+=x[j]*a[j][i];

 for(i=0;i<MAX;i++)
   f+=ftemp[i]*x[i]/2+b[i]*x[i];/*計算每次迭代時多項式的近似值*/

 for(i=0;i<MAX;i++)
   for(j=0;j<MAX;j++)
   gtemp1[i]+=a[i][j]*x[j];
 for(i=0;i<MAX;i++)
   g[i]=gtemp1[i]+b[i];     /*求多項式的梯度矩陣*/


 for(i=0;i<MAX;i++)
   for(j=0;j<MAX;j++)
   gtemp[i]+=-ani[i][j]*g[j];/*求得搜索方向*/
 
 for(i=0;i<MAX;i++)       /*求每次迭代的最優步長*/ 
   p+=g[i]*g[i];
   w=sqrt(p);             /*求梯度矩陣范式的模*/

 for(i=0;i<MAX;i++)
   for(j=0;j<MAX;j++)
   gtemp2[i]+=g[j]*a[j][i];

 for(i=0;i<MAX;i++)
   q+=gtemp2[i]*g[i];

 if(q!=0) 
 t=p/q;                     /*求得最優下降步長t=p/q*/


  /*輸出部分*/
 fprintf(fp,"第 %d 次迭代:",count);
 fprintf(fp,"目標函數值f=%lf  ",f+c);
 fprintf(fp,"梯度范式的模w=%lf  ",w);
 for(i=0;i<MAX;i++)
 fprintf(fp,"x[%d]=%lf\t",i+1,x[i]);
 fprintf(fp,"\n\n");
 /*輸出部分*/


 for(i=0;i<MAX;i++)   
 x[i]=x[i]+t*gtemp[i];      /*代入最優步長對t進行循環迭代*/

 count++;                   /*計數器*/
 
}
 
 
 fprintf(fp,"最后的結果為:\n");      /*在文件中輸出最后的結果*/
 fprintf(fp,"總共迭代%d次,控制的精度為%f.\n",count-1,ET); 
 fprintf(fp,"目標函數值f=%lf  ",f+c);
 fprintf(fp,"最優解為: ");
 for(i=0;i<MAX;i++)
 fprintf(fp,"x[%d]=%lf\t",i+1,x[i]);
 fclose(fp);

 printf("\nThe last result is:\n");  /*在TC中輸出最后的結果*/
 printf("A total of %d iterative.\n",count-1);
 printf("f=%lf  ",f+c);
 printf("The optimal solution is: ");
 for(i=0;i<MAX;i++)
 printf("x[%d]=%lf\t",i+1,x[i]);
 getch();
 
}