该楼层疑似违规已被系统折叠 隐藏此楼查看此楼

/*广义逆矩阵求解*/

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