用 c 语言实现 高斯消去法 解线性方程组解 (512 * 512 的方程组)大概需要多少时间？

[tiaoci]

在我 PIII900的机器上运行了 3-4秒，是不是太慢了？

[tiaoci]

up

[tiaoci]

竟然没人做过这种事？

[cfadongdongcfa]

差不多了，高斯的效率本来就低

[sunwebmaster]

vc中测试一段程序的运行时间我这样写了：
void CTestTimeView::OnMyTestTime()
{
CClientDC pDC(this);
CString str;
int timer=0,t,x=50,y=-10;
double counter=0.0;
SYSTEMTIME st,et;
GetSystemTime(&st);
while(counter<1.0e+6)
counter+=0.1;

GetSystemTime(&et);
if(t=et.wHour-st.wHour)
timer+=t*60*60*1000;
if(t=et.wMinute-st.wMinute)
timer+=t*60*1000;
if(t=et.wSecond-st.wSecond)
timer+=t*1000;
t=et.wMilliseconds-st.wMilliseconds;
timer+=t;
str.Format("start time=%d(h)",et.wHour);
pDC.TextOut(x,y+=20,str);
str.Format("start time=%d(h)",st.wHour);
pDC.TextOut(x,y+=20,str);
str.Format("start time=%d(m)",et.wMinute);
pDC.TextOut(x,y+=20,str);
str.Format("start time=%d(m)",st.wMinute);
pDC.TextOut(x,y+=20,str);
str.Format("start time=%d(s)",et.wSecond);
pDC.TextOut(x,y+=20,str);
str.Format("start time=%d(s)",st.wSecond);
pDC.TextOut(x,y+=20,str);
str.Format("start time=%d(ms)",et.wMilliseconds);
pDC.TextOut(x,y+=20,str);
str.Format("start time=%d(ms)",st.wMilliseconds);
pDC.TextOut(x,y+=20,str);
str.Format("total time=%d(ms)",timer);
pDC.TextOut(x,y+=20,str);


} 但是，小时显示的总比实际的少8个小时，大家看看

[寻开心]

从这个程序看, 你的代码已经写的很好了
但是也不是没有可以优化的地方

比如行交换的地方, 使用memcpy替代逐个的交换数据
例如, 在malloc 指针的时候, 再多malloc一个临时空间double * t
行交换就是
memcpy(t, i行, 长度)
memcpy(i行, j行, 长度)
memcpy(j行, t, 长度)
对max的除法转换对统一的max的倒数的乘法
在对行循环的时候, 数组的开始地址可以用指针方式

这些优化可以对效率有一定的提高, 但是不会是翻倍的提高

[寻开心]

另外, 既然你对你的算法和代码都很有信心,
那么你就应该认可你的执行速度
你在怀疑什么呢?

[DiabloWalkOnTheEarth]

已经如此简单的东西了,如果不是使用别的算法,速度都不会有明显的提高了.
swap时用 memcpy 应该没什么用,程序主要运行的时间不是在 swap 而在消去操作的地方.

[tiaoci]

你有看到C#的相应数据吗？关键是相同的代码C#也能达到3-4秒和 30秒

而且在1024*1024时，C#的速度甚至比c未优化前的还快（30秒：40秒）

如果只有这么点差异，那么用c是否还有意义？

[nicknide]

选择你相信的，相信你选择的

仅此而已。

如果你认为C#对你更加合适，大可以用C#啊……

C的优势并不在数值运算，而在于对地层的接近上，因此适合做系统。
如果搂主追求数值计算的速度，推荐用Fortain

[mathe]

增加内存，换个L2Cache大些的CPU的就可以加快速度

[dongpy]

<数值分析>里编过高斯消去法解线性方程组,不过没有512 * 512 这么大.
PIII900的机器上运行3-4秒,我觉得不是很慢.

[寻开心]

它的运算量是可以预计的啊
对于第n列来说,涉及到的是 n-1*(交换或者是消元计算, 交换也好, 消元也好都是n次计算), 所以是O(n^2)的计算两
因此总体变成上三角时候,就是 O(n^3)的计算量

再次回退计算结果的时候, 和上面的过程类似
所以总的计算量就是O(n^3)这个数量级的

如果你认为你的速度慢, 首先是检测你的计算方法有无问题
其次是看在那个部分运行占据的时间最多, 然后优化它
时间计算,如果在VC下用GetTickCount这样的毫米级的时钟就可以了

[寻开心]

单纯的对一个程序的测试估算时间没有多大的意义

重点要做的是检测算法本身以及代码编写当中可优化的地方

任何一个理论上的算法, 到具体的实践代码的转换过程当中都是有可能造成数倍,几十倍甚至成百上千倍的运行速度上的差异的, 原因就是编码水平和优化技巧的使用的问题

[tiaoci]

我的算法肯定是没有问题的，Gauss消去法是一个标准算法

要说“理论上的算法, 到具体的实践代码的转换过程当中都是有可能造成
数倍,几十倍甚至成百上千倍的运行速度上的差异的”

我可不大相信，应当不会这么夸张的，

这么悬殊的差异只有可能是算法本身造成的

要不谁来写一个解线形方程组得过程，看看是否能更快

目前在我PIII 900的机器下，使用MS的VC
计算512*512
如果不使用编译器优化大概是 3-4秒，
如果使用编译器优化 /O2大概2秒
计算1024*1024
如果不使用编译器优化大概是 40秒，
如果使用编译器优化 /O2 大概20秒

顺便贴一下使用C#实现的速度结果
计算512*512 大概是 4秒
计算1024*1024 大概是 30秒

[寻开心]

为什么不可能呢
比如在涉及到循环的算法当中
循环体内的任何一个改动, 都可以被循环本身所放大
也许这个改动只是使得算法多执行了一个乘法
那么一个100次的循环就使得这个乘法被对执行了100次
如果再有二层,三层循环, 这个执行就会被再次的放大放大

比如公式 x/y/z 和 x/(y*z) 在公式的层面上是一致的
但是在电脑的执行速度上差异是很大的
类似的问题在实现具体算法的时候,会遇到很多的
这些都是优化自己的算法实现代码的重点所在
尤其是在多重循环的最里层循环当中, 一个小的改动就可以把运行速度影响很大

[tiaoci]

to happy__888

这样的代码最多快个一倍两倍，即使某个操作被多做了100倍，

但是和 1000000( O(N^3) )的总操作量来说,几乎不算什么

而且现代的编译器已经相当智能了，象这种 x/y/z 这样的

都能自动优化掉的，现代编译器甚至能发现 /y/z 被做了多次

而使用一个预计算量来代替

[tiaoci]

看看谁能写一个解线形方程组得算法,

能比现在的2秒再快上点

[tiaoci]

这里是一个选主列的gauss过程，我没有做刻意的优化，
比方 m[i * n + j] 可以先计算 m[i * n] 然后通过指针取得
因为我相信编译器会自动优化的，看谁能改得更快点

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>

void T()
{
time_t ltime;
time(<ime);
printf("The time is %s\n", ctime( <ime ));
}

int gauss(double * m, int n, double * b)
{
int i,j,k;
double temp = 0.0;
double max = 0.0;
int pos = 0;
for(i = 0; i < n; i++)
{
max = fabs(m[i * n + i]);
pos = i;

// 选最大主列
for(j = i + 1; j < n; j++)
{
if(fabs(m[j * n +i]) > max)
{
max = fabs(m[j * n + i]);
pos = j;
}
}
if(i != pos)
{
// 交换两行的数据
for(j = i; j < n; j++)
{
temp = m[i * n + j];
m[i * n + j] = m[pos * n + j];
m[pos * n + j] = temp;
}
temp = b[i];
b[i] = b[pos];
b[pos] = temp;
}

// 当前行规一化
max = m[i * n + i];
if(fabs(max) < 0.00000001) return 0;
for(j = i + 1; j < n; j++)
{
m[i * n + j] = m[i * n + j] / max;
}
b[i] = b[i] / max;

// 执行消去操作
for(j = 0; j < n; j++)
{
if(j == i) continue;
for(k = i + 1; k < n; k++)
{
m[j * n + k] = m[j * n + k] - m[j * n + i] * m[i * n + k];
}
b[j] = b[j] - m[j * n + i] * b[i];
}
}
return 1;
}

void main()
{
int n = 1024;
double * m = (double *)malloc(sizeof(double)*n*n);
double * a = (double *)malloc(sizeof(double) * n);
int i,j;
srand(1);
for(i = 0; i < n; i++)
{
for(j = 0; j < n; j++)
{
m[i * n + j] = (double)rand() / (double)RAND_MAX;
}
a[i] = (double)rand() / (double)RAND_MAX;
}
T();
if(gauss(m, n, a) == 1)
{
T();
}
}

[gongluyang]

to tiaoci(我挑刺，我快乐)
我怎么用你这个算法得不到正确的结果呢
就算我输入a={1,1,1,-1}; b ={1,0}
应该得到的是{0.5,0.5}
怎么我得到的是{0.043620,0.078776}

还是有什么别的原因？？