我编的c语言“线性代数”处理函数，虚心接受哥哥姐姐的指教

linux_niu · 发表于 2008-7-28 11:27:20

也算是开源吧！！！！！！
刚刚学习c语言，代码难免有很多不足。
若能得到您的改进，感激不尽！！！！！
以下代码用gcc4编译都能顺利执行，并且结果正确。

代码在2楼
函数说明在3楼

linux_niu · 发表于 2008-7-28 11:28:06

#include <stdio.h>
#include <malloc.h>
#include <stdlib.h>
int cofactor(int n,double *p,int k,int l)
{
int i,j,h,m;
for(i=0,h=0;i<n;i++)
{
if(i==k)continue;
for(j=0;j<n;j++,h++)
{
if(j==l){h--;continue;}
p[h]=p[i*n+j];
}
}
m=n-1;
return m;
}
double value_det(int n,double *det)
{
int maxi,maxj,i,j;
double value=1,*p;p=NULL;
i=0;
do
{
i++;
p=(double *)malloc(n*n*sizeof(double));
if(i>100){printf("function value_det() failed because malloc()!\n");exit(0);}
}
while(p==NULL);
for(i=0;i<n*n;i++)p[i]=det[i];
while(1)
{
if(n>=3)
{
for(i=0,maxi=0,maxj=0;i<n;i++)for(j=0;j<n;j++)if(p[i*n+j]>p[maxi*n+maxj]){maxi=i;maxj=j;}
for(i=0,j=0;i<n;i++)if(p[maxi*n+i]==0)j++;
for(i=0;i<n;i++)if(p[i*n+maxj]==0)j--;
if(j>=0)
{
for(j=0;j<n;j++)
{
if(j==maxj||p[maxi*n+j]==0)continue;
for(i=0;i<n;i++)
{
if(i==maxi)continue;
p[i*n+j]=p[i*n+j]-p[i*n+maxj]*p[maxi*n+j]/p[maxi*n+maxj];
}
}
}
else {
for(i=0;i<n;i++)
{
if(i==maxi||p[i*n+maxj]==0)continue;
for(j=0;j<n;j++)
{
if(j==maxj)continue;
p[i*n+j]=p[i*n+j]-p[maxi*n+j]*p[i*n+maxj]/p[maxi*n+maxj];
}
}
}
if((maxi+maxj)%2==1)value*=-1;
value*=p[maxi*n+maxj];
n=cofactor(n,p,maxi,maxj);
}
else if(n==2){value*=(p[0]*p[3]-p[1]*p[2]);break;}
else if(n==1){value*=p[0];break;}
}
free(p);
return value;
}
double *solve_cramer(int n,double *mtr)
{
double *x,*p,d;
int i=0,j,k;
do
{i++;
x=(double *)malloc(n*sizeof(double));
if(i>100){printf("function solve_cramer() failed because malloc()!\n");exit(0);}
}
while(x==NULL);i=0;
do
{i++;
p=(double *)malloc(n*n*sizeof(double));
if(i>100){printf("function solve_cramer() failed because malloc()!\n");exit(0);}
}
while(p==NULL);
for(j=0,k=0;k<n*n;j++){if((j+1)%(n+1)==0)continue;p[k]=mtr[j];k++;}
d=value_det(n,p);
if(d==0){free(x);x=NULL;}
else
{
for(i=0;i<n;i++)
{
for(j=0,k=0;k<n*n;j++){if((j+1)%(n+1)==0)continue;p[k]=mtr[j];k++;}
for(k=0;k<n;k++)p[k*n+i]=mtr[k*(n+1)+n];
x[i]=value_det(n,p)/d;
}
}
free(p);
return x;
}
double *multiply_mtr(double *mtr1,double *mtr2,int m,int s,int n)
{
double *p1,*p2,*p;
int i=0,k,j;
do
{
i++;
p1=(double *)malloc(m*s*sizeof(double));
if(i>100){printf("function multiply_mtr() failed because malloc()!\n");exit(0);}
}
while(p1==NULL);i=0;
do
{
i++;
p2=(double *)malloc(n*s*sizeof(double));
if(i>100){printf("function multiply_mtr() failed because malloc()!\n");exit(0);}
}
while(p2==NULL);i=0;
do
{
i++;
p=(double *)malloc(m*n*sizeof(double));
if(i>100){printf("function multiply_mtr() failed because malloc()!\n");exit(0);}
}
while(p==NULL);
for(i=0;i<m*s;i++)p1[i]=mtr1[i];
for(i=0;i<s*n;i++)p2[i]=mtr2[i];
for(i=0;i<m*n;i++)p[i]=0;
for(i=0;i<m;i++)
{
for(j=0;j<n;j++)
{
for(k=0;k<s;k++)p[i*n+j]+=p1[i*s+k]*p2[k*n+j];
}
}
free(p1);
free(p2);
return p;
}
double *transp_mtr(double *mtr,int m,int n)
{
double *p;
int i=0,j;
do
{i++;
p=(double *)malloc(m*n*sizeof(double));
if(i>100){printf("function transp_mtr() failed by malloc()!\n");exit(0);}
}
while(p==NULL);
for(i=0;i<n;i++)
for(j=0;j<m;j++)
p[i*m+j]=mtr[j*n+i];
return p;
}
double *inv_mtr(double *mtr,int n)
{
int i=0,j,k;
double *inver,*p,a;
inver=NULL;
p=NULL;
do
{i++;
inver=(double *)malloc(n*n*sizeof(double));
if(i>100){printf("function inv_mtr() failed by malloc()!\n");exit(0);}
}
while(inver==NULL);i=0;
do
{i++;
p=(double *)malloc(n*n*sizeof(double));
if(i>100){printf("function inv_mtr() failed by malloc()!\n");exit(0);}
}
while(p==NULL);
for(i=0;i<n*n;i++)p[i]=mtr[i];
if(value_det(n,p)==0){printf("no inverse of the matrix\n");free(p);free(inver);inver=NULL;}
else
{
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
for(k=0;k<n*n;k++)p[k]=mtr[k];
k=cofactor(n,p,i,j);
inver[i*n+j]=value_det(k,p);
if((i+j)%2==1)inver[i*n+j]*=-1;
}
}
free(p);
p=inver;
inver=transp_mtr(p,n,n);
free(p);
a=value_det(n,mtr);
for(i=0;i<n*n;i++)inver[i]/=a;
}
return inver;
}

复制代码

linux_niu · 发表于 2008-7-28 11:28:48

int cofactor(int n,double *p,int k,int l)
求n阶行列式除去k行l列後的余子式（处理时注意下标）。传送数组的首地址给double *p。返回处理後行列式的阶数。
===========================================================
double value_det(int n,doulbe *p)
求n阶行列式的值。传送数组的首地址给double *p。返回行列式的值。
===========================================================
double *solve_cramer(int n,double *mtr)
用克拉默法则解含有n个未知数的线性方程组，即系数行列式为n阶。将存放增广矩阵的数组首地址传送给double *mtr，返回存放n个解的数组的首地址。若无解返回NULL。
===========================================================
double *multiply_mtr(double *mtr1,double *mtr2,int m,int s,int n)
计算m*s与s*n两矩阵相乘，将存放两矩阵的数组的首地址分别传给double *mtr1和double *mtr2。返回计算後的矩阵首地址。
===========================================================
double *transp_mtr(double *mtr,int m,int n)
矩阵转置。将存放m*n矩阵的数组的首地址传送给double *mtr。返回存放转置後的矩阵的首地址。
===========================================================
double *inv_mtr(double *mtr,int n)
求n阶矩阵的逆矩阵，返回存放逆矩阵的数组首地址。若无逆矩阵返回NULL。
n>=2

remote fish · 发表于 2008-7-28 12:13:41

提个代码风格的建议, 个人认为 do-while 循环写成楼主这样有点别扭. 很多人喜欢写成下面的形式

do {
} while ();

复制代码

重点在于, while() 和 } 写在一行让人一眼就看得出它不是一个 while 循环, 而是 do-while 循环

remote fish · 发表于 2008-7-28 12:20:02

另外, 在 malloc 失败的时候进行循环重试似乎意义并不是太大. 根据 man malloc 来看, malloc 并不会受到 signal 干扰而失败, 唯一提到的可能的失败情况是 ENOMEM, 这种情况下, 重试也没什么用

remote fish · 发表于 2008-7-28 12:24:24

像这样的批量赋值

for(i=0;i<m*s;i++)p1[i]=mtr1[i];

复制代码

如果用块拷贝的方式会有更高的效率, 毕竟是数学库嘛, 效率还是比较重要的. 或者说, 有什么别的需要考虑的因素吗?

memcpy(p1, mtr1, m * s * sizeof(*mtr1));

复制代码

linux_niu · 发表于 2008-7-28 12:34:47

Post by remote fish;1879410
像这样的批量赋值
for(i=0;i<m*s;i++)p1[i]=mtr1[i];
复制代码
如果用块拷贝的方式会有更高的效率, 毕竟是数学库嘛, 效率还是比较重要的. 或者说, 有什么别的需要考虑的因素吗?
memcpy(p1, mtr1, m * s * sizeof(*mtr1));
复制代码

这个太好了，
我都不知道有这个函数，去学习下！！！
十分感谢

这个就是开源的好处，呵呵

remote fish · 发表于 2008-7-28 12:40:00

回得真快. 楼主再看看 man memset 的说明吧

linux_niu · 发表于 2008-7-28 12:45:55

Post by remote fish;1879415
回得真快. 楼主再看看 man memset 的说明吧

是在shell下执行$man memset吗？
就是这样了！！！

xin@debian:~$ man memset
No manual entry for memset
xin@debian:~$

remote fish · 发表于 2008-7-28 13:15:42

debian 默认没有装 glibc 的 manpages
手工装一下吧, 好像是 manpages-dev 这个包. 记不清了

sudo aptitude install manpages-dev

复制代码

		自动登录	找回密码
密码			注册