又一些c语言“线性代数”函数(有更新,并有一些疑问向大家请教,永不过期)

linux_niu · 发表于 2008-8-1 10:21:27

这个函数的数组循环什么的有什么毛病吗？
还是别的地方有毛病？？
为什么有的矩阵计算正确，而有的计算错误？
double **rdechl_mtr(double **mtr,int m,int n)

请到二楼看代码，
三楼有函数说明
谢谢您了！！！！

linux_niu · 发表于 2008-8-1 10:22:08

#include <stdio.h>
#include <malloc.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
double **make_matrix(int m,int n)
{
double **p;
int i,k,h;
k=0;
do
{
p=(double **)malloc(m*sizeof(double *));
k++;
if(k>100){printf("Function make_matrix() failed!\n");exit(0);}
}while(p==NULL);
k=0;
for(i=0;i<m;i++)p[i]=NULL;
do
{
for(i=0,h=0;i<m;i++)
{
if(p[i]==NULL)p[i]=(double *)malloc(n*sizeof(double));
if(p[i]==NULL)h++;
}
k++;
if(k>100){printf("Function make_matrix() failed!\n");exit(0);}
}while(h!=0);
return p;
}
double **make_parr(int m)
{
double **p;
int k;
k=0;
do
{
p=(double **)malloc(m*sizeof(double *));
k++;
if(k>100){printf("Funtion make_parr() failed!\n");exit(0);}
}while(p==NULL);
return p;
}
double *make_array(int n)
{
double *p;
int k;
k=0;
do
{
p=(double *)malloc(n*sizeof(double));
k++;
if(k>100){printf("Function make_array() failed!\n");exit(0);}
}while(p==NULL);
return p;
}
void free_matrix(double **mtr,int m)
{
int i;
for(i=0;i<m;i++)free(mtr[i]);
free(mtr);
}
void val_cpy(double **mtr1,double **mtr2,int m,int n)
{
int i;
for(i=0;i<m;i++)memcpy(mtr1[i],mtr2[i],n*sizeof(double));
}
void idt_mtr(double **mtr,int n)
{
int i,j;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
if(i==j)mtr[i][j]=1;
else mtr[i][j]=0;
}
}
void exch_mtr(double **mtr,int n,int r,int s)
{
double *t;
t=make_array(n);
memcpy(t,mtr[r],n*sizeof(double));
memcpy(mtr[r],mtr[s],n*sizeof(double));
memcpy(mtr[s],t,n*sizeof(double));
}

复制代码

有问题函数:

double **rdechl_mtr(double **mtr,int m,int n)
{
double **p,tmp;
int i,j,k,max,zero,r;
p=make_matrix(m,n);
val_cpy(p,mtr,m,n);
for(j=0,k=0;j<n;j++,k++)
{
for(i=k,max=k,zero=k;i<m;i++)
{
if(fabs(p[i][j])>fabs(p[max][j]))max=i;
if(p[i][j]==0)zero++;
}
if(zero==m){k--;continue;}
else
{
if(max!=k)exch_mtr(p,n,k,max);
for(i=j,tmp=p[k][j];i<n;i++)p[k][i]/=tmp;
for(i=0;i<m;i++)
{
if(i==k)continue;
for(r=j,tmp=p[i][j];r<n;r++)p[i][r]-=p[k][r]*tmp;
}
}
}
return p;
}

复制代码

linux_niu · 发表于 2008-8-1 10:22:49

doulbe **make_matrix(int m,int n)
动态创建m*n二维数组。返回数组首地址。
===========================================================
double **make_parr(int m)
动态创建存放m个指针的指针数组。
===========================================================
double *make_array(int n)
动态创建有n个元素的一维数组。返回数组首地址。
===========================================================
void free_matrix(double **mtr,int m)
释放由malloc()创建的有m行的二维数组
===========================================================
void val_cpy(double **mtr1,double **mtr2,int m,int n)
将首地址为double **mtr2的m*n数组复制到首地址为double **mtr1的的数组中。（只要数组循环不越界）
===========================================================
void idt_mtr(double **mtr,int n)
将n阶矩阵初始化成单位阵
===========================================================
void exch_mtr(double **mtr,int n,int r,int s)
交换有n列矩阵的第r行与第s行
===========================================================
double **rdechl_mtr(double **mtr,int m,int n)
求矩阵mtr的行最简型。
返回计算後的矩阵首地址

asert · 发表于 2008-8-3 23:56:53

看楼主学线性代数学的那么好，以后有线性代数方面的问题我PM你哦

remote fish · 发表于 2008-8-4 10:51:37

楼主指的运算错误是说结果不对还是程序异常? 给出一组可以导致错误的数据吧, 这让别人才好测试

linux_niu · 发表于 2008-8-5 06:29:43

利用rdechl_mtr()计算矩阵的行最简型,
结果正确的是:

m=4
n=5
Input 4*5 Matrix:
2 -1 -1 1 2
1 1 -2 1 4
4 -6 2 -2 4
3 6 -9 7 9
the 4*5 Matrix:
2.0 -1.0 -1.0 1.0 2.0
1.0 1.0 -2.0 1.0 4.0
4.0 -6.0 2.0 -2.0 4.0
3.0 6.0 -9.0 7.0 9.0
reduced echelon matrix:
1.0 0.0 -1.0 0.0 4.0
0.0 1.0 -1.0 0.0 3.0
0.0 0.0 0.0 1.0 -3.0
0.0 0.0 0.0 0.0 0.0

复制代码

结果错误的是

Matrix m*n:
m=4
n=5
Input 4*5 Matrix:
1 -1 3 -4 3
3 -3 5 -4 1
2 -2 3 -2 0
3 -3 4 -2 -1
the 4*5 Matrix:
1.0 -1.0 3.0 -4.0 3.0
3.0 -3.0 5.0 -4.0 1.0
2.0 -2.0 3.0 -2.0 0.0
3.0 -3.0 4.0 -2.0 -1.0
reduced echelon matrix:
1.0 -1.0 0.0 0.0 0.0
0.0 0.0 1.0 0.0 0.0
0.0 0.0 0.0 1.0 0.0
0.0 0.0 0.0 0.0 1.0

复制代码

上面代码的正确结果应该是

1.0 -1.0 0.0 2.0 -3.0
0.0 0.0 1.0 -2.0 2.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0

复制代码

linux_niu · 发表于 2008-8-5 06:31:04

Post by remote fish;1882300
楼主指的运算错误是说结果不对还是程序异常? 给出一组可以导致错误的数据吧, 这让别人才好测试

fish要是能帮我找出原因,真感激不尽啊!!!都看一天多了也没找到原因,就放弃了

remote fish · 发表于 2008-8-5 10:49:24

用 gdb 跟着跑了一遍, 其实错误很明显的, 是浮点数精度的问题.

当 rdechl_mtr 函数中 j 达到 3 的时候, p 中保存的已经是结果式了, 然后在对 0 进行计数时, linux_niu 兄是用

if (p[i][j] == 0) zero++;

复制代码

这样的语句. 如果 p[j] 是整形量, 那么这样完全正确, 但是不幸的是, p[j] 是浮点型, 而且其结果又是经过几轮相减后得到的一个接近 0 的小量, 但是又不是 0, 这就是浮点数误差.

在 j == 3 的时候, p[2][j] 的值是 -5.921e-16, 实际上就是零了, 但是又无法满足 p[2][j] == 0 的判断

这种情况下, 我们可以把此语句改成

if (fabs(p[i][j]) < 0.0001) zero++;

复制代码

0.0001 也可以换成一个稍大一点或稍小一点的值, 只要保证 [1 它远小于可能出现在矩阵中的最小的非零绝对值], [2 又远大于可能达到的累积误差], 这两个条件即可.

实际上, 以楼主的目的来看, 用分数来表示一个数可能更合理一些, 不过这要进行一些先期的封装.

P.S. 建议 linux_niu 兄对代码进行以下几点改进:
* 对关键变量的作用进行说明
* 在关键代码处加上注释
* 用空行分隔语义段
* for(r=j,tmp=p[j];r<n;r++)p[r]-=p[k][r]*tmp; 这样的代码写成

tmp = p[i][j];

for (r = j; r < n; r++)

p[i][r] -= p[k][r] * tmp;

复制代码

的形式更易于阅读
* 统一函数接口风格. 比如将名字都改成 matrix_make, matrix_free, matrix_exch, matrix_rdechl, array_make, 这样一看就知道哪些是 matrix 相关的, 哪些是 array 相关的
* malloc 分配的内存要记得释放, 比如 exch_mtr() 中的 t
* 更根本的, 在一个函数中申请一块内存给外部使用, 比如 rdechl_mtr() 中的 p, 这种方式非常不好, 说到底还是代码结构的问题. 在另一帖中 realtang 等兄弟们已经推荐过一些好的资料了, 楼主不妨参考一下.

linux_niu · 发表于 2008-8-6 17:07:34

十分感谢9楼
您的留言我都认真学习了，
谢谢！！！！！

LiYanrui · 发表于 2008-8-6 20:28:41

如果真得需要许多线性代数运算，自己又觉得实在太累，可以考虑用 gsl 吧。如果要学习，也可以看看 gsl 源码。

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