C_and_CPP 版 (精华区)

发信人: hittih (hittih), 信区: C_and_CPP
标  题: Re: 请问哪里有矩阵求逆的C  函数，谢谢！
发信站: 哈工大紫丁香 (Fri Mar  5 10:58:44 2004), 站内信件

         int matrix_athwart(double a[],int n)
{ 
int *is,*js,i,j,k,l,u,v;
double d,p;
is=new int[n];
js=new int[n];
for (k=0; k<=n-1; k++)
{ d=0.0;
for (i=k; i<=n-1; i++)
for (j=k; j<=n-1; j++)
{ l=i*n+j; p=fabs (a[l]);
if (p>d) { d=p; is[k]=i; js[k]=j;}
}
if (d+1.0==1.0)

{ delete[] is; is=NULL; delete[] js; js=NULL; printf("err**not inv\n");

return(0);
}
if (is[k]!=k)
for (j=0; j<=n-1; j++)
{ u=k*n+j; v=is[k]*n+j;
p=a[u]; a[u]=a[v]; a[v]=p;
}
if (js[k]!=k)
for (i=0; i<=n-1; i++)
{ u=i*n+k; v=i*n+js[k];
p=a[u]; a[u]=a[v]; a[v]=p;
}
l=k*n+k;
a[l]=1.0/a[l];
for (j=0; j<=n-1; j++)
if (j!=k)
{ u=k*n+j; a[u]=a[u]*a[l];}
for (i=0; i<=n-1; i++)
if (i!=k)
for (j=0; j<=n-1; j++)
if (j!=k)
{ u=i*n+j;
a[u]=a[u]-a[i*n+k]*a[k*n+j];
}
for (i=0; i<=n-1; i++)
if (i!=k)
{ u=i*n+k; a[u]=-a[u]*a[l];}
}
for (k=n-1; k>=0; k--)
{ if (js[k]!=k)
for (j=0; j<=n-1; j++)
{ u=k*n+j; v=js[k]*n+j;
p=a[u]; a[u]=a[v]; a[v]=p;
}
if (is[k]!=k)
for (i=0; i<=n-1; i++)
{ u=i*n+k; v=i*n+is[k];
p=a[u]; a[u]=a[v]; a[v]=p;
}
}
delete[] is;  is=NULL;
delete[] js; js=NULL;
return(1);
}
【 在 shyfox 的大作中提到: 】
: 请高手指教，谢谢！

--

※ 来源:．哈工大紫丁香 bbs.hit.edu.cn [FROM: 210.46.69.10]