该楼层疑似违规已被系统折叠 隐藏此楼查看此楼
/*广义逆矩阵求解*/
double *ginv(double *A,int m,int n)
{
double *An,*R,*R1,*R2,*R3,*C,*C1,*C3,*B,b,d,temp,*X,*Y,*Z,*Bn;
int i,j,bn,k,x,y;
An=(double *)malloc(n*m*(sizeof(double)));
R=(double *)malloc(m*m*(sizeof(double)));
C=(double *)malloc(n*n*(sizeof(double)));
B=(double *)malloc(n*m*sizeof(double));
Bn=(double *)malloc(n*m*sizeof(double));
R1=(double *)malloc(m*m*sizeof(double));
R2=(double *)malloc(m*m*sizeof(double));
R3=(double *)malloc(m*m*sizeof(double));
C1=(double *)malloc(n*n*sizeof(double));
C3=(double *)malloc(n*n*sizeof(double));
X=(double *)malloc(max(n,m)*max(n,m)*sizeof(double));
Y=(double *)malloc(max(n,m)*max(n,m)*sizeof(double));
Z=(double *)malloc(max(n,m)*max(n,m)*sizeof(double));
//定义单位阵R
for(x=0;x
{
for(k=0;k
{
if(x==k)
{
R[x*m+k]=1;
}
else
{
R[x*m+k]=0;
}
}
}
//定义单位阵C
for(x=0;x
{
for(k=0;k
{
if(k==x)
{
C[x*n+k]=1;
}
else
{
C[x*n+k]=0;
}
}
}
//定义矩阵b=a
for(i=0;i
{
for(j=0;j
B[i*n+j]=A[i*n+j];
}
i=0;
bn=n;
while(i
{
//初始化
b=B[i*n+i];
j=i;
for(k=i;k
{
if(fabs(B[k*n+i])>fabs(b))
{
b=B[k*n+i];
j=k;
}
}
if(b==0&&i==bn-1)
{
break;
}
else if(b==0&&i!=bn-1)
{
for(k=0;k
{
temp=B[k*n+i];
B[k*n+i]=B[k*n+bn-1];
B[k*n+bn-1]=temp;
}
for(x=0;x
{
for(k=0;k
{
if(x==k)
C1[x*n+k]=1;
else
C1[x*n+k]=0;
}
}
C1[i*n+(bn-1)]=1;
C1[(bn-1)*n+i]=1;
C1[i*n+i]=0;
C1[(bn-1)*n+(bn-1)]=0;
mult(C,C1,X,n,n,n);
for(x=0;x
{
for(k=0;k
{
C[x*n+k]=X[x*n+k];
}
}
bn--;
}
else if(j!=i)
{
for(k=0;k
{
temp=B[j*n+k];
B[j*n+k]=B[i*n+k];
B[i*n+k]=temp;
}
for(x=0;x
{
for(y=0;y
{
if(x==y)
R1[x*m+y]=1;
else
R1[x*m+y]=0;
}
}
R1[i*m+j]=1;
R1[j*m+i]=1;
R1[i*m+i]=0;
R1[j*m+j]=0;
mult(R1,R,X,m,m,m);
for(x=0;x
{
for(k=0;k
{
R[x*m+k]=X[x*m+k];
}
}
}
if(b!=1&&b!=0)
{
for(k=0;k
B[i*n+k]=B[i*n+k]/b;
}
for(x=0;x
{
for(y=0;y
{
if(x==y)
R2[x*m+y]=1;
else