一个小型矩阵库

原创 2006年05月19日 11:10:00

                                                                        一个小型矩阵库

本文是一个小型的矩阵库,对做一些数学应用开发很有用 !!

#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

FILE *fp;
/****************************************************************************************
*
*                         关于一些用到的数学函数
*
****************************************************************************************/

/***************************************************************************************
*              
*                            数学函数--求整数的幂
*
****************************************************************************************/
int ipow(int i,int j)
{
 int k,s=1;
 for(k=1;k<=j;k++)
  s=s*i;
 return s;
}
/***************************************************************************************
*              
*                            数学函数--数组
*
****************************************************************************************/
int Kronecker(int x,int y,int z)
{
 return (x*y==z)?1:0;
}

/***************************************************************************************
*                              <ds_matrix.h文件>
*                           关于矩阵及其操作函数
*
*        double Get(matrix m,int i,int j) ;          求矩阵特定位置元素值       
*        int Set(pmatrix pm,int i,int j,const double d);       给矩阵特定位置元素负值
*        void DisplayM(matrix m);                   输出矩阵
*        matrix Sum(matrix m1,matrix m2);           求矩阵和
*        matrix ZhuanZhi(matrix m);       矩阵转置
*        double Det(matrix m);              用全选主元消去法求行列式的值
*        matrix ShuCheng(matrix m,double d);       实现矩阵数乘
*        matrix Multi(matrix m_A,matrix m_B);           实现矩阵乘法
*        matrix ZLJ(const matrix m1,const matrix m2);      矩阵的张量积
*        matrix Inv(matrix from);                        矩阵求逆
*        double NeiJi(matrix m1,matrix  m2);    求两个向量的内积
*        void FreeM(pmatrix pm)  ;               释放矩阵内存        
****************************************************************************************/

/***************************************************************************************
*              
*                           定义数据结构---矩阵
*
****************************************************************************************/

typedef struct tag_matrix
{
  int hang;
  int lie;
 double * p;
}matrix;
typedef matrix* pmatrix;

/*声明用到的变量和函数*/
extern FILE *fp;
extern matrix I(int N);
/***************************************************************************************
*              
*                            求矩阵特定位置元素值
*
****************************************************************************************/


double Get(matrix m,int i,int j)
{
 if(m.p==NULL||i>=m.hang||j>=m.lie)

 {
  printf("不能从矩阵中正确取值,可能是所给参数有误 !!/n");
        fprintf(fp,"不能从矩阵中正确取值,可能是所给参数有误 !!/n");
  return 0;
 }
 return *(m.p+i*m.lie+j);
}
/***************************************************************************************
*              
*                            给矩阵特定位置元素负值
*
****************************************************************************************/
int Set(matrix m,int i,int j,const double d)
{
 if(m.p==NULL||i>=m.hang||j>=m.lie)
 {
  printf("不能正确付值,可能是所给参数有误 !!/n");
  fprintf(fp,"不能正确付值,可能是所给参数有误 !!/n");
  return 0;
 }
   *(m.p+i*m.lie+j)=d;
 return 1;

}
/***************************************************************************************
*              
*                            输出矩阵
*
****************************************************************************************/
void DisplayM(matrix m)
{
 int i,j;
    fprintf(fp,"/n输出矩阵%d行%d列:/n",m.hang,m.lie);
 for(i=0;i<m.hang;i++)
 {

  for(j=0;j<m.lie;j++)
  { 
            fprintf(fp," %e",Get(m,i,j));
  }
        fprintf(fp,";/n");
 }

}
void FreeM(pmatrix pm)
{
 if(pm->p!=NULL)
 {
  pm->hang=0;
  pm->lie=0;
  free(pm->p);
  pm->p=NULL;
 }
}

/***************************************************************************************
*              
*                             求矩阵和
*
****************************************************************************************/
matrix Sum(matrix m1,matrix m2)
{
 int i,j;
 double d;
    matrix m;
    m.p=NULL;
 if(m1.hang!=m2.hang||m1.lie!=m2.lie)
 {
  printf("所给矩阵行列数不相等,不能求和!!");
  fprintf(fp,"所给矩阵行列数不相等,不能求和!!");
  return m;
 }
 m.hang=m1.hang;
 m.lie=m1.lie;
 m.p=(double *)malloc(sizeof(double)*m.hang*m.lie);
 for(i=0;i<m.hang;i++)
  for(j=0;j<m.lie;j++)
  {

   d=Get(m1,i,j)+Get(m2,i,j);
   Set(m,i,j,d);
  }
     return m;
}

/***************************************************************************************
*              
*                           矩阵转置
*
****************************************************************************************/
matrix ZhuanZhi(matrix m)
{
 int i,j;
 matrix tm;
    tm.hang=m.lie;
 tm.lie=m.hang;
 tm.p=(double*)malloc(sizeof(double)*tm.hang*tm.lie);
 for(i=0;i<tm.hang;i++)
  for(j=0;j<tm.lie;j++)
  {
   Set(tm,i,j,Get(m,j,i));
  }
    return tm;

}

/***************************************************************************************
*              
*                            用选列元消去法求行列式的值
*
****************************************************************************************/
double Det(matrix m)

 int i,j,row,col,k;
 double max,temp,det=1,n=m.hang;
 matrix m1;
    if(m.hang!=m.lie)
 {
  printf("所给矩阵非方阵,不能计算行列式值/n");
  fprintf(fp,"所给矩阵非方阵,不能计算行列式值/n");
  return 0;
 }
     m1.hang=m.hang;
  m1.lie=m.lie;
  m1.p=(double*)malloc(sizeof(double)*m1.hang*m1.lie);
  for(i=0;i<m1.hang;i++)
   for(j=0;j<m1.lie;j++)
   {
    Set(m1,i,j,Get(m,i,j));
   }
 for(k=0;k<n;k++)
 {
  /*找主元*/
  max=0;row=col=k;
  for(i=k;i<n;i++)
   for(j=k;j<n;j++)
   {
    temp=fabs(Get(m1,i,j));
    if(max<temp){max=temp;row=i;col=j;}

   }
   /*交换行列,将主元交换到k行,k列*/
   if(max==0)return 0;
   if(row!=k)
   {
    for(j=0;j<n;j++)
    {
     temp=*(m1.p+row*m1.lie+j);
     Set(m1,row,j,Get(m1,k,j));
                    Set(m1,k,j,temp);
    }
   }
   if(col!=k)
   {
    for(i=0;i<n;i++)
    {
     temp=Get(m1,i,k);
     Set(m1,i,k,Get(m1,i,col));
                    Set(m1,i,col,temp);
    }
   }
   det*=Get(m1,k,k);
   for(j=k+1;j<n;j++)Set(m1,k,j,Get(m1,k,j)*1.0/Get(m1,k,k));
   for(j=k+1;j<n;j++)
    for(i=k+1;i<n;i++)
     *(m1.p+i*m1.lie+j)-=Get(m1,i,k)*Get(m1,k,j);
   for(i=0;i<k;i++)Set(m1,i,k,0);
            Set(m1,k,k,1);
 }
 FreeM(&m1);
   return det;

}
/***************************************************************************************
*              
*                                    实现矩阵数乘
*
****************************************************************************************/
matrix ShuCheng(matrix m,double d)
{
 int i,j;
 matrix tm;
 tm.p=(double*)malloc(sizeof(double)*m.hang*m.lie);
 if(tm.p==NULL)
 {
  printf("申请内存失败/n");
  fprintf(fp,"申请内存失败/n");
  return tm;
 }
 tm.hang=m.hang;
 tm.lie=m.lie;
 for(i=0;i<m.hang;i++)
  for(j=0;j<m.lie;j++)
   Set(tm,i,j,d*Get(m,i,j));
 return tm;
}
/***************************************************************************************
*              
*                                    实现矩阵内积
*
****************************************************************************************/
double NeiJi(matrix m1,matrix  m2)
 {
  int i;
  double s=0;
  if(m1.hang==1)/*横向量*/
  {
   if(m2.hang!=1)
   {
    printf("/n所给向量不能进行内积运算!!/n");
    fprintf(fp,"/n所给向量不能进行内积运算!!/n");
    return 0;
   }
     for( i=0;i<m1.lie;i++)
          s+=Get(m1,0,i)*Get(m2,0,i);
   return s;
  }
  else if(m1.lie==1)/*纵向量*/
  {
  if(m2.lie!=1)
   {
    printf("/n所给向量不能进行内积运算!!/n");
    fprintf(fp,"/n所给向量不能进行内积运算!!/n");
    return 0;
   }
  for( i=0;i<m1.hang;i++)
          s+=Get(m1,i,0)*Get(m2,i,0);
  return s;
  }
  else
  {
             printf("/n所给向量不能进行内积运算!!/n");
    fprintf(fp,"/n所给向量不能进行内积运算!!/n");
    return 0;
  }


 }
/***************************************************************************************
*              
*                                实现向量的范数
*
****************************************************************************************/
double Norm(matrix m)
{

 double d;
 d=NeiJi(m,m);
 d=sqrt(d);
 return d;
}
matrix Sub(matrix m1,matrix m2)
{
 int i,j;
 double d;
 matrix m;
 m.p=NULL;
 if(m1.hang!=m2.hang||m1.lie!=m2.lie)
 {
  printf("/n所给矩阵不能进行减法运算/n");
  fprintf(fp,"/n所给矩阵不能进行减法运算/n");
  return m;
 }
 m.hang=m1.hang;
 m.lie=m1.lie;
    m.p=(double *)malloc(sizeof(double)*m.hang*m.lie);
 for(i=0;i<m.hang;i++)
  for(j=0;j<m.lie;j++)
  {

   d=Get(m1,i,j)-Get(m2,i,j);
   Set(m,i,j,d);
  }
    return m;
}
/***************************************************************************************
*              
*                                    实现矩阵乘法
*
****************************************************************************************/
matrix Multi(matrix m_A,matrix m_B)
{   
 int i,j,k;
 double t;
 matrix m_C;
 m_C.p=NULL;
 if(m_A.lie!=m_B.hang)
 {
  printf("您给的矩阵不能进行乘法运算/n");
  fprintf(fp,"您给的矩阵不能进行乘法运算/n");
  return m_C;
 }

 m_C.hang=m_A.hang;
 m_C.lie=m_B.lie;
 m_C.p=(double *)malloc(sizeof(double)*m_C.hang*m_C.lie);
 if(m_C.p==NULL)
 {
  printf("申请内存失败/n");
  fprintf(fp,"申请内存失败/n");
  return m_C;
 }
 for(i=0;i<m_C.hang;i++)   /*计算C的i+1行j+1列元素*/
  for(j=0;j<m_C.lie;j++)
  {
   t=0;
            for(k=0;k<m_A.lie ;k++)
   {
    t+=Get(m_A,i,k)*Get(m_B,k,j);
                Set(m_C,i,j,t);
   }
  }

    return m_C;
}
/***************************************************************************************
*              
*                          求矩阵的张量积
*
****************************************************************************************/
matrix ZLJ( matrix m1, matrix m2)
{
 int i,j;
 double d;
 matrix m;
 m.hang=m1.hang*m2.hang;
 m.lie=m1.lie*m2.lie;
 m.p=(double *)malloc (sizeof(double)*m.hang*m.lie); 
 if(m.p==NULL)
 {
  printf("申请内存失败/n");
  fprintf(fp,"申请内存失败/n");
  return m;
 }
    for(i=0;i<m.hang;i++)
  for(j=0;j<m.lie;j++)
  {
   d=Get(m1,i/m2.hang,j/m2.lie)*Get(m2,i%m2.hang,j%m2.lie);
            Set(m,i,j,d);
  }
 return m;
}
/***************************************************************************************
*              
*                                   矩阵求逆
*
****************************************************************************************/
matrix Inv(matrix from_m)

 int i,j,row,col,k,n,*p,t2,t1;
 double max,temp;
 matrix to,from;
    to.p=NULL;
    if(Det(from_m)==0)return  to; /*此矩阵无逆*/
 from=I(from_m.hang);
 for(i=0;i<from.hang;i++)
  for(j=0;j<from.lie;j++)
   Set(from,i,j,Get(from_m,i,j));
 to=I(from.hang);
    n=from.hang;
    p=(int *)malloc(sizeof(int)*n);
 for(i=0;i<n;i++)
 p[i]=i;
 for(k=0;k<n;k++)
 {
  /*找主元*/
  max=0;
  row=k;
  col=k;
   for(i=k;i<n;i++)
 
   {
    temp=fabs(Get(from,i,k));
    if(max<temp)
    {
     max=temp;
              row=i;
    
    }
   }
   /* 交换行列,将主元调整到k行k列上*/
   if(row!=k)
   {
    for(j=0;j<n;j++)
    {
     temp=Get(from,row,j);
                    Set(from,row,j,Get(from,k,j));
                    Set(from,k,j,temp);

                    temp=Get(to,row,j);
                    Set(to,row,j,Get(to,k,j));
                    Set(to,k,j,temp);
    }

    for(i=0;i<n;i++)
    {
     if(k==p[i])
      t1=i;
     if(row==p[i])
      t2=i;
    }
   
   }
  
   /*k行归一*/
   
   for(j=k+1;j<n;j++)
   {
    temp=Get(from,k,j)*1.0;
    temp=temp/Get(from,k,k);
    Set(from,k,j,temp);
   }
   for(j=0;j<n;j++)  
   {
    temp=Get(to,k,j)*1.0;
    temp=temp/Get(from,k,k);
    Set(to,k,j,temp);
   }
            Set(from,k,k,1);

   for(j=k+1;j<n;j++)
   {
    for(i=0;i<k;i++)
    {
     temp=Get(from,i,j)-Get(from,i,k)*Get(from,k,j);
     Set(from,i,j,temp);
    }
                 for(i=k+1;i<n;i++)
     {
                    temp=Get(from,i,j)-Get(from,i,k)*Get(from,k,j);
     Set(from,i,j,temp);
     }
   }

    for(j=0;j<n;j++)
    {
     for(i=0;i<k;i++)
     {
                        temp=Get(to,i,j)-Get(from,i,k)*Get(to,k,j);
         Set(to,i,j,temp);
     }
     for(i=k+1;i<n;i++)
     {
                        temp=Get(to,i,j)-Get(from,i,k)*Get(to,k,j);
         Set(to,i,j,temp);
     }
    }

   for(i=0;i<n;i++)
    Set(from,i,k,0);
   Set(from,k,k,1);
 }
 FreeM(&from);
    return to;
}


double MaxT(matrix m)/*采用幂法求最大特征值*/
{
 int i,k;
 double m0,m1,t;
 matrix x0,x1;
 if(m.hang!=m.lie)
 {
  printf("/n所给矩阵不能求条件数!!/n");
  fprintf(fp,"/n所给矩阵不能求条件数!!/n");
  return 0;
 }
 x0.hang=m.hang;
 x0.lie=1;
 x0.p=(double*)malloc(sizeof(double)*x0.hang);
 for(i=0;i<x0.hang;i++)  /*给x0赋初值*/
 {
  Set(x0,i,0,0.5);
 }
 k=0;                  /*给m0赋值*/
 for(i=1;i<x0.hang;i++)
 {
  if(fabs(Get(x0,i,0))>fabs(Get(x0,k,0)))
   k=i;
 }
 m0=Get(x0,k,0);
   for(i=1;i<x0.hang;i++)
 {
  Set(x0,i,0,(1.0*Get(x0,i,0))/m0);
 }
 while(1)
 {
  x1=Multi(m,x0);
        k=0;                 
     for(i=1;i<x1.hang;i++)
  {
     if(fabs(Get(x1,i,0))>fabs(Get(x1,k,0)))
      k=i;
  }
     m1=Get(x1,k,0);
     for(i=0;i<x1.hang;i++)
  {
   t=(Get(x1,i,0)*1.0)/m1;
            Set(x0,i,0,t);
  }
     if(fabs(m1-m0)<1e-6)
  {
      FreeM(&x1);
      FreeM(&x0);
         return m1;
  }
     m0=m1;
 }

}
double Cond(matrix m) /*求条件数*/
{
 double t1,t2;
 matrix tm1,tm2,tm3;
 tm1=ZhuanZhi(m);
 tm2=Multi(tm1,m);
 t1=MaxT(tm2);
 FreeM(&tm1);
 FreeM(&tm2);
    tm1=Inv(m);
 tm2=ZhuanZhi(tm1);
 tm3=Multi(tm2,tm1);
 t2=MaxT(tm3);
 FreeM(&tm1);
    FreeM(&tm2);
 FreeM(&tm3);
 return sqrt(t1*t2);
}

/***************************************************************************************
*              
*                          <construct_matrix.h文件>
*                             构造一些用到的矩阵
*
*
*              matrix Zeros(int hang,int lie);     构造全零矩阵
*              matrix I(int n) ;        构造单位矩阵
*              matrix Trig(double x,double y,double z,int N); 构造三对角矩阵
*              matrix X(int i,int j,int sig,int tau,int ep,int m,int N) ; 构造大X矩阵
*              matrix Xpiao(int i,int j,int sig,int tau,int N); 构造X一弯矩阵
*         matrix  S(int N);   构造S矩阵
*              matrix DeltaXX(int N); 构造DeltaXX矩阵
*              matrix DeltaYY(int N); 构造DeltaYY矩阵
*              matrix DeltaZZ(int N); 构造DeltaZZ矩阵
*
*
*
****************************************************************************************/


/***************************************************************************************
*              
*                             构造全零矩阵
*
****************************************************************************************/
matrix Zeros(int hang,int lie)
{
      int i,j;
   matrix m;
      m.p=(double *)malloc(sizeof(double)*hang*lie);
      if(m.p==NULL)
   {
     printf("/n申请内存失败/n");
        return m;
   }
   m.hang=hang;
   m.lie=lie;
   for(i=0;i<hang;i++)
    for(j=0;j<lie;j++)
      Set(m,i,j,0);
   return m;
}
/***************************************************************************************
*              
*                             构造单位矩阵
*
****************************************************************************************/
matrix I(int n)          
{
 int i,j;
 matrix m;
    m.p=(double *)malloc(sizeof(double)*n*n);
    if(m.p==NULL)
 {
     printf("/n申请内存失败/n");
        return m;
  }
  m.hang=n;
  m.lie=n;
   for(i=0;i<n;i++)
    for(j=0;j<n;j++)
    {
     if(i==j)
             Set(m,i,j,1);
     else
                     Set(m,i,j,0);
    }
    return m;
}
/***************************************************************************************
*              
*                         构造三对角矩阵
*
****************************************************************************************/
matrix Trig(double x,double y,double z,int N)
{
 int i,j;
 matrix m;
 m=Zeros(N,N);
 for(i=1;i<N-1;i++)
  for(j=1;j<N-1;j++)
  { if(i==j)
   {
    Set(m,i,j,y);
       Set(m,i,j-1,x);
    Set(m,i,j+1,z);
   }
  }
 Set(m,0,0,y);
 Set(m,0,1,z);
 Set(m,N-1,N-2,x);
 Set(m,N-1,N-1,y);
 m.hang=N;
 m.lie=N;
 return m;
}
//使用BI——CG方法求解方程,方程为:AX=B
matrix BI_CG(matrix A,matrix B,int N,int *Num,double* t)
{
 FILE *tfp1,*tfp2,*tfp3;
 int i,len;
 double ZhiShu1,ZhiShu2,tempZhiShu;
 matrix x1,x2,p1,p2,pW1,pW2,r1,r2,rW1,rW2,v,tm1,tm2,U,tempx,Cham,tS;
 double alpha,beta,rho1=0,rho2,start,finish;

    tfp1=fopen("time.txt","w");
    tfp2=fopen("Zhishu.txt","w");
    tfp3=fopen("b.txt","w");


 start=clock();
 len=(2*N-1)*(2*N-1)*(2*N-1);
 x1.hang=(2*N-1)*(2*N-1)*(2*N-1);
 x1.lie=1;
 x1.p=(double *)malloc (sizeof(double)*x1.hang*x1.lie);
    p1.hang=(2*N-1)*(2*N-1)*(2*N-1);
 p1.lie=1;
 p1.p=(double *)malloc (sizeof(double)*p1.hang*p1.lie);
    p2.hang=(2*N-1)*(2*N-1)*(2*N-1);
 p2.lie=1;
 p2.p=(double *)malloc (sizeof(double)*p2.hang*p2.lie);
 pW1.hang=(2*N-1)*(2*N-1)*(2*N-1);
 pW1.lie=1;
 pW1.p=(double *)malloc (sizeof(double)*pW1.hang*pW1.lie);
 pW2.hang=(2*N-1)*(2*N-1)*(2*N-1);
 pW2.lie=1;
 pW2.p=(double *)malloc (sizeof(double)*pW2.hang*pW2.lie);
    r1.hang=(2*N-1)*(2*N-1)*(2*N-1);
 r1.lie=1;
 r1.p=(double *)malloc (sizeof(double)*r1.hang*r1.lie);
 r2.hang=(2*N-1)*(2*N-1)*(2*N-1);
 r2.lie=1;
 r2.p=(double *)malloc (sizeof(double)*r2.hang*r2.lie);
    rW1.hang=(2*N-1)*(2*N-1)*(2*N-1);
 rW1.lie=1;
 rW1.p=(double *)malloc (sizeof(double)*rW1.hang*rW1.lie);
 rW2.hang=(2*N-1)*(2*N-1)*(2*N-1);
 rW2.lie=1;
 rW2.p=(double *)malloc (sizeof(double)*rW2.hang*rW2.lie);
    v.hang=(2*N-1)*(2*N-1)*(2*N-1);
 v.lie=1;
 v.p=(double *)malloc (sizeof(double)*v.hang*v.lie);
                             /*赋初值*/
 for(i=0;i<x1.hang;i++)
  Set(x1,i,0,0.5);
    tm1=Multi(A,x1);
 tm2=Sub(B,tm1);
 FreeM(&tm1);
 for(i=0;i<p1.hang;i++)
 {
  Set(p1,i,0,Get(tm2,i,0));
 }
 for(i=0;i<r1.hang;i++)
 {
  Set(r1,i,0,Get(tm2,i,0));
 } 
 for(i=0;i<rW1.hang;i++)
 {
  Set(rW1,i,0,Get(tm2,i,0));
 }
    for(i=0;i<pW1.hang;i++)
 {
  Set(pW1,i,0,Get(tm2,i,0));
 }
 FreeM(&tm2);
    rho1=NeiJi(rW1,r1);
 *Num=0;
 x2.p=NULL;


 U=left_u(N);
    tS=S(N);

ZhiShu2=10;
    while(1)
 {

  tm1=Multi(A,p1);/*1--v*/
        for(i=0;i<len;i++)   
  {
   Set(v,i,0,Get(tm1,i,0));
  }
  FreeM(&tm1);
        alpha=rho1*1.0/(NeiJi(v,pW1)); /*2---alpha*/
  tm1=ShuCheng(p1,alpha);   /*3---x*/
  FreeM(&x2);
  x2=Sum(x1,tm1);
  FreeM(&tm1);
  tm1=Sub(x2,x1);
  (*Num)++;
  //      if(Norm(tm1)<1e-6)
//   break;

     tempx=Multi(tS,x2);//begin
  Cham=Sub(tempx,U);
  FreeM(&tempx);
  ZhiShu1=log10(Norm(Cham));
  fprintf(tfp2,"%f;/n",ZhiShu1);
     finish=clock();
     *t=(finish-start)*1.0/CLK_TCK;
        fprintf(tfp1,"%f;/n",*t);

     if(ZhiShu1<-12.0)
   break;
  tempZhiShu=fabs(ZhiShu1-ZhiShu2);
  if(tempZhiShu<1e-6)
  {
   fprintf(tfp2,"误差趋于稳定/n");
   break;
  }
      ZhiShu2=ZhiShu1;


                                             //end


  FreeM(&tm1);
  for(i=0;i<len;i++)
   Set(x1,i,0,Get(x2,i,0));

  tm1=ShuCheng(v,alpha);  /*4---r*/
  tm2=Sub(r1,tm1);
  FreeM(&tm1);
        for(i=0;i<len;i++)
   Set(r2,i,0,Get(tm2,i,0));
  FreeM(&tm2);
  for(i=0;i<len;i++)
   Set(r1,i,0,Get(r2,i,0));
  tm1=ZhuanZhi(A);  /*5---rW*/
        tm2=ShuCheng(tm1,alpha );
  FreeM(&tm1);
  tm1=Multi(tm2,pW1);
  FreeM(&tm2);
  tm2=Sub(rW1,tm1);
  FreeM(&tm1);
        for(i=0;i<len;i++)
          Set(rW2,i,0,Get(tm2,i,0));
  for(i=0;i<len;i++)
   Set(rW1,i,0,Get(rW2,i,0));
  FreeM(&tm2);
  rho2=NeiJi(r2,rW2);  /*6---rho*/
  beta=rho2*1.0/rho1;/*7---beta*/
  rho1=rho2;
  tm1=ShuCheng(p1,beta);/*8---p*/
  tm2=Sum(r2,tm1);
  FreeM(&tm1);
  for(i=0;i<len;i++)  
   Set(p2,i,0,Get(tm2,i,0));
  for(i=0;i<len;i++)
   Set(p1,i,0,Get(p2,i,0));
  FreeM(&tm2);
        tm1=ShuCheng(pW1,beta); /*9---pW*/
  tm2=Sum(rW2,tm1);
  FreeM(&tm1);
  for(i=0;i<len;i++)
   Set(pW2,i,0,Get(tm2,i,0));
         for(i=0;i<len;i++)
   Set(pW1,i,0,Get(tm2,i,0));
  FreeM(&tm2);
}
 FreeM(&v);
 FreeM(&x1);
 FreeM(&r1);
    FreeM(&r2);
 FreeM(&rW1);
 FreeM(&rW2);
 FreeM(&p1);
 FreeM(&p2);
 FreeM(&pW1);
 FreeM(&pW2);
 finish=clock();
 *t=(finish-start)*1.0/CLK_TCK;
 return x2;

}
/*使用共轭梯度法求线形方程组的解*/
matrix CG(matrix A,matrix B,int N,int *Num,double* t)
{
    int i;
 matrix x0,x1,p0,p1,r0,r1,tm1,tm2;
 double alpha ,beta;
 time_t start,finish;
 start=clock();
 x0.hang=(2*N-1)*(2*N-1)*(2*N-1);
 x0.lie=1;
 x0.p=(double *)malloc (sizeof(double)*x0.hang*x0.lie);
    p0.hang=(2*N-1)*(2*N-1)*(2*N-1);
 p0.lie=1;
 p0.p=(double *)malloc (sizeof(double)*p0.hang*p0.lie);
 p1.hang=(2*N-1)*(2*N-1)*(2*N-1);
 p1.lie=1;
 p1.p=(double *)malloc (sizeof(double)*p1.hang*p1.lie);
    r0.hang=(2*N-1)*(2*N-1)*(2*N-1);
 r0.lie=1;
 r0.p=(double *)malloc (sizeof(double)*r0.hang*r0.lie);
 r1.hang=(2*N-1)*(2*N-1)*(2*N-1);
 r1.lie=1;
 r1.p=(double *)malloc (sizeof(double)*r1.hang*r1.lie);

 /*赋初值*/

 for(i=0;i<x0.hang;i++)
  Set(x0,i,0,0.5);
 tm1=Multi(A,x0);
 tm2=Sum(B,tm1);
 FreeM(&tm1);
 for(i=0;i<x0.hang;i++)
  Set(x0,i,0,Get(tm2,i,0));
     for(i=0;i<r0.hang;i++)
  Set(r0,i,0,Get(tm2,i,0));
  FreeM(&tm2);
  *Num=0;
  /*迭代*/
  x1.p=NULL;
  while(1)
  {
   tm1=Multi(A,p0);
   alpha=(NeiJi(p0,r0)*1.0)/(NeiJi(tm1,p0));
   FreeM(&tm1);
   tm1=ShuCheng(p0,alpha);
   FreeM(&x1);
   x1=Sum(x0,tm1);
   FreeM(&tm1);
   (*Num)++;
   tm1=Sub(x1,x0);
   if(Norm(tm1)<1e-6)
   {
    FreeM(&tm1);
    break;
   }
   FreeM(&tm1);
         for(i=0;i<x1.hang;i++)
    Set(x0,i,0,Get(x1,i,0));
          tm1=Multi(A,x1);
    tm2=Sub(B,tm1);
    FreeM(&tm1);
    for(i=0;i<r0.hang;i++)
    {
    Set(r1,i,0,Get(tm2,i,0));
    Set(r0,i,0,Get(tm2,i,0));
    }
       FreeM(&tm2);
    tm1=Multi(A,r1);
    tm2=Multi(A,p0);
    beta=(-1.0*NeiJi(tm1,p0))/(NeiJi(tm2,p0));
    FreeM(&tm1);
    FreeM(&tm2);
    tm1=ShuCheng(p0,beta);
    tm2=Sum(r1,tm1);
    FreeM(&tm1);
       for(i=0;i<p1.hang;i++)
  {
   Set(p0,i,0,Get(tm2,i,0));
   Set(p1,i,0,Get(tm2,i,0));
  }
    FreeM(&tm2);
  }
  FreeM(&tm1);
  FreeM(&tm2);
  FreeM(&r0);
  FreeM(&r1);
     FreeM(&p0);
     FreeM(&p1);
  FreeM(&x0);
  finish=clock();
  *t=(finish-start)*1.0/CLK_TCK;
 
  return x1;
}
matrix NI(matrix A,matrix B)/*使用直接求逆法求线形方程组的解*/
{
 matrix m,tm;
 tm=Inv(A);
 m=Multi(tm,B);
 FreeM(&tm);
 return m;
}

版权声明:本文为博主原创文章,未经博主允许不得转载。

相关文章推荐

网络编程(编写一个小型服务器)--csapp11章

网络应用随处可见,表现的形式也各不相同,有趣的是,在不同的表现形式之下都是基于相同的编程模型,依赖于相同的编程接口,因此学习网络还是比较保值的,因为这么多基础设备还在运行着,基本机制也在短时间内很难改...

如何使用AsyncSocket编写一个非常小型的群聊系统

本人刚刚从事ios开发半年, 前几天想了解一下苹果的Socket编程, 然后去百度东西, 发现许多人都是建议使用AsyncSocket, 然后就去了解了一下, 现在我把我的学习经验总结一下, 哪里说的...

打造自己的专属Linux(一):快速建立一个小型Linux

跟着马哥学Linux已经快一个月了,在前几天创建一个小型Linux,大小不足10M,萌生了一个想法:DIY一个自己的Linux,所以准备写一个系列篇,通过宿主机的方式,一步一步创建一个小型Linux,...

一个小型字符串工具类MString

自己实现了一个字符串工具类,不怎么完善,凑合这用, 里面的韩硕或者变量都还比较简单,很容易理解 #ifndef MSTRING_H #define MSTRING_H /********...

如何建一个小型网站

你是一个技术人员,你想创建一个站点,或许是一时心血来潮,或许是为了自己的兴趣爱好,或许是...「创业」前的热身准备?那么,如何少走弯路的构建起来你的 Web 站点呢?别笑,不是每个做技术的都捣鼓过个人...

C#开发一个小型编译器

设计任务为开发一个小型编译器,主要包括4个部件:词法分析器、语法分析器、语义分析器、目标代码生成器。 一、课程设计的目的与要求 1.通过设计、编写和调试词法分析程序,了解词法扫描器的组成结构、不同种...

综合图形学算法及OpenGL技术开发一个小型图形软件系统

开发一个图形软件包,具有画直线、曲线、曲面填充、几何变换、裁剪等功能:

大二的时候的一个小项目---实现四则运算的逻辑分析实现(小型编译器)

package com.my.calculator; public class Number { private double number = 0.0; public Number(doub...

初学如何构建一个小型项目的流程与心得

这篇博客我来说说我的看法,对于怎样合理的

一个小型调试器的原理与设计

(按:由于鄙人文字奇烂,字迹尤丑,有不通凝塞之处还请看官忍耐。小文一篇,聊引方家解颐) 由于之前一直对gdb的原理感到一定兴趣,所以最近小研究了一下调试器相关的东西(Linux平台)。gdb有一篇...
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:深度学习:神经网络中的前向传播和反向传播算法推导
举报原因:
原因补充:

(最多只允许输入30个字)