一个小型矩阵库
来源：互联网 发布：js random int 编辑：程序博客网 时间：2024/06/11 16:53
一个小型矩阵库
本文是一个小型的矩阵库，对做一些数学应用开发很有用 ！！
#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;
}
- 一个小型矩阵库
- Dillo-一个小型网页浏览器
- 一个小型的搜索引擎
- 一个小型的操作系统
- 实现一个小型通讯录
- 一个小型名字字典
- 一个小型课程表
- 小型库
- Thinking in C++读书笔记--2.2一个小型的C库
- sChart.js：一个小型简单的图表库
- 一个小型的企业信息项目
- 一个小型的溢出实验
- 如何建一个小型网站
- C#开发一个小型编译器
- 一个小型的安装脚本
- Cage 是一个 Java 实现的验证码图片生成库，快速、小型和简单。
- 一个小型VC项目的开发
- 一个小型管理系统的pb实现
- JBuilder 2005开发Applet游戏全接触
- 垃圾回收
- DOM
- linux文件目录结构介绍
- dataReader转化为dataTable
- 一个小型矩阵库
- 修改SourceInsight3.5试用版
- ASP中添加记录并返回ID的方法
- Meta标签详解，在网上转的，希望对大家有用
- 怎么实现在FireFox IE Opera Safari 都可以正常播放WMV和MOV的网页播放器代码
- SqlServer存储过程编写经验和优化
- Spring--简单使用quartz实现定时作业
- 使用Spring邮件抽象层发送邮件
- 7 Reasons Why Web Apps Fail
