• 设为首页
  • 点击收藏
  • 手机版
    手机扫一扫访问
    迪恩网络手机版
  • 关注官方公众号
    微信扫一扫关注
    公众号
登陆 注册

【转】vs2010+matlab求解线性方程组

原作者: [db:作者] 来自: [db:来源] 收藏 邀请

代码转帖：http://www.cnblogs.com/youthlion/archive/2012/05/29/2524599.html

最近仿08年TOG上一篇骨架提取的文章Skeleton extraction by mesh contraction，其中涉及到线性方程组的最小二乘解问题，即Ax = b。

最开始使用了Armadillo库进行求解，程序写完后发现矩阵A的规模与顶点数的平方成正比，不使用稀疏矩阵的话只能计算很小的模型，但Armadillo没有提供稀疏矩阵模块。听说Eigen库有稀疏矩阵模块，又查了下Eigen库，但是发现Eigen库的稀疏矩阵求解线性方程组的功能只能用于A为方阵的情况。最后考虑用matlab结合vs2010的方式实现。

Armadillo和Eigen都是很绿色的线性代数库，都是以泛型编程的方式实现。Armadillo非常简练，文档也小巧精悍，上手很快，底层依赖lapack和blas库，环境配置方式写在了之前的一篇文档里http://www.cnblogs.com/youthlion/archive/2012/05/15/2501465.html。Eigen库更加重量级一些，功能更加全面，文档详细，不依赖于任何其他底层库。环境配置也很简单，不多说了。下面主要记录一下vs2010和matlab混用的方法。

其实matlab提供了多种工具，既可以在matlab中调用其他语言写好的模块，也可以在其他语言中调用matlab生成的模块。因为不熟悉matlab，所以我选择的方式是用matlab生成动态链接库，在c/c++中调用。这是Matlab Compiler提供的功能。

vs2010环境配置方式如下：

1、首先要在vs2010中设置所需头文件的路径，在我的电脑上是D:\Program Files\MATLAB\R2011a\extern\include；

2、指定库文件路径，D:\Program Files\MATLAB\R2011a\extern\lib\win32\microsoft；

3、在linker中添加附加依赖库mclmcrrt.lib，Over。

然后在matlab中写好需要的功能，第一次用matlab，代码蠢笨，大家见笑了。

首先，要在matlab中实现稀疏矩阵的创建（文件名createsparse.m），其中m，n分别是行数和列数：

function A = createsparse(m, n)
A = sparse(m,n);
end

接下来要实现对稀疏矩阵中的非零元素赋值（文件名eleassign.m），其中i，j分别是行索引和列索引，value是元素的值：

function A = eleassign( A, i, j, value )
A(i,j)=value;
end

矩阵元素赋值完成后就是线性方程组的求解。这个函数是直接在mathworks官网上下载的，代码比较长，这里只给出声明（文件名lsmr.m）：

function [x]...
   = lsmr(A, b, lambda, atol, btol, conlim, itnlim)

实际上原始的输入和输出参数要更多一些，由于我只用到一个输出，即方程组的解向量，所以把其他输出删掉了。

现在有了这三个函数，分别对应三个.m文件，下面需要用这三个函数生成对应的动态链接库，在matlab中使用如下命令：

mcc -W cpplib:SparseSolver -T link:lib createsparse.m eleassign.m lsmr.m

其中SparseSolver是库文件名，-T link:lib 的作用是指定输出文件类型为库文件。

matlab busy一会儿以后，就会产生一系列文件，用到的是其中的三个：SparseSolver.h、SparseSolver.lib和SparseSolver.dll，把这三个文件拷贝到工程文件夹中。

其中SparseSolver.lib是导入库，在vs2010的工程属性中的附加依赖库中添加这个文件。

SparseSolver.h是生成的c/c++函数声明，比如上面实现的几个函数，在头文件中分别为：

extern LIB_SparseSolver_CPP_API void MW_CALL_CONV createsparse(int nargout, 
                                                                mwArray& A, 
                                                                const mwArray& m, 
                                                                const mwArray& n);
extern LIB_SparseSolver_CPP_API void MW_CALL_CONV eleassign(int nargout, 
                                                            mwArray& A, 
                                                            const mwArray& A_in1, 
                                                            const mwArray& i, 
                                                            const mwArray& j, 
                                                            const mwArray& value);
extern LIB_SparseSolver_CPP_API void MW_CALL_CONV lsmr(int nargout, 
                                                        mwArray& x, 
                                                        const mwArray& A, 
                                                        const mwArray& b, 
                                                        const mwArray& lambda, 
                                                        const mwArray& atol, 
                                                        const mwArray& btol, 
                                                        const mwArray& conlim, 
                                                        const mwArray& itnlim);　

在这些函数中，第一个参数是输出参数的个数。其他都是实际的输入输出参数。C++中与matlab互相传数据都是用这种mwArray的数据类型。比如要C++中需要向matlab中传递一个向量，可能要这样：

...
double vec[] = {1,2,3};
mwArray _vec(3,1,mxDOUBLE_CLASS,mxREAL);
_vec.SetData(vec,3);
...

mwArray _vec(3,1,...) 中，3是行数，1是列数。

有两件事要注意，一是matlab中矩阵的存储是列优先存储，比如矩阵

在matlab中存储的顺序为1，1，1，1，-1，1。而在C++中我们习惯的存储方式是{1,1,1,-1,1,1}，所以在用SetData拷贝数据的时候必须先转换行列顺序。

二是向matlab传数据只能是以mwArray传，即使只是一个普通数值，也必以mwArray的方式传过去，比如这样：

...
double somevalue[] = {1};
mwArray _somevalue(1,1,mxDOUBLE_CLASS,mxREAL);
_somevalue.SetData(somevalue,1);
...

当然mwArray也支持随机访问，比如：

...
mwArray _somevalue(1,1,mxDOUBLE_CLASS,mxREAL);
_somevalue(1,1) = 5.0f;
...

确实不太方便，不过比起寻找各种库所花费的力气，值了。

下面是一些功能测试代码，很乱：

#include <iostream>
#include "SparseSolver.h"
 
using namespace std;
int main()
{
    //mclmcrInitialize();
    if (!mclInitializeApplication(NULL,0)) 
    {
        std::cerr << "could not initialize the application properly"
            << std::endl;
        return -1;
    }
    if( !SparseSolverInitialize() )
    {
        std::cerr << "could not initialize the library properly"
            << std::endl;
        return -1;
    }
    else
    {
        try
        {
            //parameters initialization
            double lamda[] = {0.0f};
            double atol[] = {0.0f};
            double btol[] = {0.0f};
            double conlim[] = {0.0f};
            double itnlim[] = {50.0f};
            mwArray _lamda(1,1,mxDOUBLE_CLASS,mxREAL);
            //_lamda.SetData(lamda,1);
            _lamda(1,1) = lamda[0];
            mwArray _atol(1,1,mxDOUBLE_CLASS,mxREAL);
            //_atol.SetData(atol,1);
            _atol(1,1) = atol[0];
            mwArray _btol(1,1,mxDOUBLE_CLASS,mxREAL);
            //_btol.SetData(btol,1);
            _btol(1,1) = btol[0];
            mwArray _conlim(1,1,mxDOUBLE_CLASS,mxREAL);
            //_conlim.SetData(conlim,1);
            _conlim(1,1) = conlim[0];
            mwArray _itnlim(1,1,mxDOUBLE_CLASS,mxREAL);
            //_itnlim.SetData(itnlim,1);
            _itnlim(1,1) = itnlim[0];
            //data
            double matA[] = {1.0,1.0,1.0,-1.0,1.0,1.0};
            double vecb[] = {2.0,1.0,3.0};
            mwArray _vecb(3,1,mxDOUBLE_CLASS,mxREAL);
            _vecb.SetData(vecb,3);
            
            //create matrix A
            mwArray _matA;
            double rows[] = {3.0f};
            double cols[] = {2.0f};
            mwArray _rows(1,1,mxDOUBLE_CLASS,mxREAL);
            mwArray _cols(1,1,mxDOUBLE_CLASS,mxREAL);
            _rows.SetData(rows,1);
            _cols.SetData(cols,1);
            createsparse(1,_matA,_rows,_cols);
            cout<<_matA<<endl;
            for(int i = 0; i < 3; i++)
            {
                for(int j = 0; j<2; j++)
                {
                    double temp[1];
                    temp[0] = matA[i*2+j];
                    mwArray _temp(1,1,mxDOUBLE_CLASS,mxREAL);
                    _temp.SetData(temp,1);
                    double indexr[1],indexc[1];
                    indexr[0] = i+1;
                    indexc[0] = j+1;
//                     mwArray _indexr(1,1,mxDOUBLE_CLASS,mxREAL);
//                     mwArray _indexc(1,1,mxDOUBLE_CLASS,mxREAL);
//                     _indexr.SetData(indexr,1);
//                     _indexc.SetData(indexc,1);
//                     eleassign(1,_matA,_matA,_indexr,_indexc,_temp);
                    _matA(i+1,j+1) = matA[i*2+j];
                }
            }
            cout<<_matA<<endl;
            int cnt = 0;
            for(int i = 0; i<3;i++)
            {
                for(int j= 0;j<2;j++)
                    _matA(i+1,j+1) = cnt++;
            }
            cout<<_matA<<endl;
            for(int i=0;i<3;i++)
            {
                for(int j=0;j<2;j++)
                {
                    double temp = _matA(i+1,j+1);
                    if(temp == 0)
                    {
                        cout<<"zero pos:"<<endl;
                        cout<<i+1<<" "<<j+1<<endl;
                    }
                }
            }
            mwArray _matB;
            double rowsb[] = {4.0f};
            double colsb[] = {4.0f};
            mwArray _rowsb(1,1,mxDOUBLE_CLASS,mxREAL);
            mwArray _colsb(1,1,mxDOUBLE_CLASS,mxREAL);
            _rowsb.SetData(rowsb,1);
            _colsb.SetData(colsb,1);
            createsparse(1,_matB,_rowsb,_colsb);
            double rbegin[] = {1.0f};
            double rend[] = {3.0f};
            double lbegin[] = {1.0f};
            double lend[] = {2.0f};
            mwArray _rbegin(1,1,mxDOUBLE_CLASS,mxREAL);
            mwArray _rend(1,1,mxDOUBLE_CLASS,mxREAL);
            mwArray _lbegin(1,1,mxDOUBLE_CLASS,mxREAL);
            mwArray _lend(1,1,mxDOUBLE_CLASS,mxREAL);
            _rbegin.SetData(rbegin,1);
            _rend.SetData(rend,1);
            _lbegin.SetData(lbegin,1);
            _lend.SetData(lend,1);
            cout<<"B"<<endl;
            cout<<_matB<<endl;
            submat(1,_matB,_matB,_rbegin,_rend,_lbegin,_lend,_matA);
            cout<<"A:"<<endl;
            cout<<_matA<<endl;
            cout<<"B"<<endl;
            cout<<_matB<<endl;
            double temp[1];
            double ok = 2;
            temp[0] = ok;
            mwArray _temp(1,1,mxDOUBLE_CLASS,mxREAL);
            _temp.SetData(temp,1);
            matmult(1,_matB,_matB,_temp);
            cout<<"B after mult"<<endl;
            cout<<_matB<<endl;
            cout<<_matB(4,4)<<endl;
            _matB(4,4) = 5.0;
            cout<<_matB(4,4)<<endl;
            mwArray _x;
            lsmr(1,_x,_matA,_vecb,_lamda,_atol,_btol,_conlim,_itnlim);
            cout<<_x<<endl;
        }
        catch (const mwException& e)
        {
            cerr << e.what() << endl;
            system("pause");
            return -2;
        }
        catch (...)
        {
            cerr << "Unexpected error thrown" << endl;
            system("pause");
            return -3;
        }     
        // Call the application and library termination routine
        SparseSolverTerminate();
    }
    mclTerminateApplication();
    system("pause");
    return 0;
}

最后总结一下，说matlab运算速度慢，效率低的，可能要么是大牛程序员，要么是人云亦云。低水平的coder，像我这种，自己写的数值分析算法肯定不如matlab实现的。所以，对于一般程序员来说，用matlab做一个后台的计算工具是一件很爽的事情。另外matlab本身有成熟的交流平台，file exchange有大量的代码资源可以下载，也能节约不少搜集资源的时间。

就写到这里，关于matlab计算效率的问题，等把骨架提取的代码改完后再补充吧。

鲜花
握手
雷人
路过
鸡蛋
该文章已有0人参与评论

请发表评论

全部评论

专题导读
More+
上一篇：
delphiXEBerlinReadProcessMemoryWriteProcessMemory发布时间：2022-07-18
下一篇：
RxLib 2.75 for Delphi D5-DX10.4 Sydney发布时间：2022-07-18
热门推荐
More+
热门话题
More+
阅读排行榜

扫描微信二维码

查看手机版网站

随时了解更新最新资讯

139-2527-9053

在线客服（服务时间 9:00～18:00）

在线QQ客服
地址：深圳市南山区西丽大学城创智工业园
电邮：jeky_zhao#qq.com
移动电话：139-2527-9053

Powered by 互联科技 X3.4© 2001-2213 极客世界.|Sitemap