本帖最后由 kerbcurb 于 2010-12-21 13:15 编辑
阵求解器
0）
用于求解M s = b这样的线性方程组，row表示方阵M的阶数。数据类型双精度。矩阵可以ROW X ROW存储即多行多列存储，也可以1 X ROW^2也就是1行多列存储。
void Solvef ( float * * M, int row, float *b, float * & s );
//------------------------------------------------------------------------------------------------------------------------------
void Solved ( double * * M, int row,double * b, double * & s );
//------------------------------------------------------------------------------------------------------------------------------
void Solveld ( long double * * M, int row,long double * b, double * & s );
//------------------------------------------------------------------------------------------------------------------------------
void Solvedl ( double * M, int row,double * b, double * & s );
//----------------------------------------------------------------------------------------------------------------------------
1）
工程文件需要包含：#include "windows.h"，然后把MatrixSolve.dll文件复制到你的工程文件夹下。在需要使用求解器的地方写上：
typedef void ( *TSolve) ( double*, int, double*, double *& );
HINSTANCE handle;
FARPROC lpFarProc;
TSolve Solve;
handle= LoadLibrary( "MatrixSolve.dll");
if( handle == NULL )
{
FreeLibrary ( handle );
}
else
{
lpFarProc = GetProcAddress ( handle, "Solvedl" );
//使用Solvedl还是Solved或者其他的，看你的矩阵存储类型
Solve = TSolve ( lpFarProc );
// A X = B，方阵A的阶数是ROW，X是待求未知数
Solve( A, ROW, B, X );
FreeLibrary ( handle );
}
2）
如果生成文件时提示fatal error LNK1107:…….，
是FAT32系统分区的问题，如果你的硬盘式用NTFS就没有问题。。。解决方法：项目à项目属性à配置属性à清单工具à使用FAT32解决方法为 是
3）
我一般使用多行多列存储系数矩阵，这里的例子使用1行多列存储，似乎二者并不影响求解时间。我感觉多行多列存储系数矩阵稍微快一点，但是有的人的说法与我恰好相反，不知道他们使用的是什么求解器。
本帖最后由 kerbcurb 于 2010-12-20 11:09 编辑
我的空间里面的：
http://www.mwtee.com/home-space-uid-18294-do-blog-id-1203.html 介绍过另一种求解器，比这个慢70倍左右。以下是他的代码
- //--------------------------------------------------------------
- int Gauss ( double a[], double b[], int n )
- {
- int *js, l, k, i, j, is, p, q;
- double d, t;
- js = ( int* ) malloc ( n * sizeof ( int ) );
- l = 1;
- for ( k = 0; k <= n - 2; k++ )
- {
- d = 0.0;
- for ( i = k; i <= n - 1; i++ )
- for ( j = k; j <= n - 1; j++ )
- {
- t = fabs ( a[i*n+j] );
- if ( t > d )
- {
- d = t;
- js[k] = j;
- is = i;
- }
- }
- if ( d + 1.0 == 1.0 ) l = 0;
- else
- {
- if ( js[k] != k )
- for ( i = 0; i <= n - 1; i++ )
- {
- p = i * n + k;
- q = i * n + js[k];
- t = a[p];
- a[p] = a[q];
- a[q] = t;
- }
- if ( is != k )
- {
- for ( j = k; j <= n - 1; j++ )
- {
- p = k * n + j;
- q = is * n + j;
- t = a[p];
- a[p] = a[q];
- a[q] = t;
- }
- t = b[k];
- b[k] = b[is];
- b[is] = t;
- }
- }
- if ( l == 0 )
- {
- free ( js );
- printf ( "Gauss funtion failed 1...\n" );
- return ( 0 );
- }
- d = a[k*n+k];
- for ( j = k + 1; j <= n - 1; j++ )
- {
- p = k * n + j;
- a[p] = a[p] / d;
- }
- b[k] = b[k] / d;
- for ( i = k + 1; i <= n - 1; i++ )
- {
- for ( j = k + 1; j <= n - 1; j++ )
- {
- p = i * n + j;
- if ( a[i*n+k] != 0 && a[k*n+j] != 0 )
- a[p] -= a[i*n+k] * a[k*n+j];
- }
- if ( a[i*n+k] != 0 && b[k] != 0 )
- b[i] -= a[i*n+k] * b[k];
- }
- }
- d = a[ ( n-1 ) *n+n-1];
- if ( fabs ( d ) + 1.0 == 1.0 )
- {
- free ( js );
- printf ( "Gauss funtion failed 2...\n" );
- return ( 0 );
- }
- b[n-1] = b[n-1] / d;
- for ( i = n - 2; i >= 0; i-- )
- {
- t = 0.0;
- for ( j = i + 1; j <= n - 1; j++ )
- if ( a[i*n+j] != 0 && b[j] != 0 )
- t += a[i*n+j] * b[j];
- b[i] = b[i] - t;
- }
- js[n-1] = n - 1;
- for ( k = n - 1; k >= 0; k-- )
- if ( js[k] != k )
- {
- t = b[k];
- b[k] = b[js[k]];
- b[js[k]] = t;
- }
- free ( js );
- return ( 1 );
- }
本帖最后由 kerbcurb 于 2010-12-20 23:46 编辑
回复 cem-uestc 的帖子
与那个代码有很大区别，他的那个代码在本主题帖子上面也贴出来了，他的代码应该有很多地方可以优化
测试过的，请给出参考意见，以及与其他求解器的比较结果等。目前还没做好写一个稀疏矩阵存储以及求解的方案，筹划中...
回复 kerbcurb 的帖子
呵呵，我们做计算电磁学的人，都是每人负责一块的，就是有人专门研究求解器的，所以效率高点，其实国外有很多开源库的，已经做得很好的，如果好好找找就会发现……
回复 kerbcurb 的帖子
呵呵，这个不是偷偷的问题，这个就是术业有专攻的问题。人家就是研究这个的东西的，况且这个是开源，不存在版权问题。如果你什么东西都自己编，那你的代码效率肯定不会很高，毕竟你要把各部分都做好是不可能的事情。
如果什么都用别人的，估计什么也做不好！现在好的东西太多了，未见几个在这些好的基础上做出什么更好的成绩，国内计算电磁学界有多少不是在重复外人做过的东西呢，课题重复也没有什么，是不是比原先的更先进？
发布一款矩阵求解器: MatrixSolve.dll
