谈谈利用开源代码来编写求解器

caoer

20# bbssbb

可以考虑考虑petsc,

northwindgod

同意楼主观点，国内搞有限元软件开发的本身就基本上是游击队，人手有限，再撒胡椒面儿把精力花到求解器这一块儿实在没太大必要。
另外别的方面也应该有侧重点，走差异化路线才有前途。

jdbear

读了楼主的帖子，我不得不说话，simwea是个开放的论坛，又见有想法、有经历、有水平的帖子，希望搞基础力学的多多聆听和学习！

pasuka

谢谢12楼提供的介绍稀疏矩阵存贮格式的文章。

英文的那篇感觉和现有的研究工作有一些重复，可能它们也要独立开发软件包吧，但愿不要低水平重复。中文的那篇感觉意义一般。

我建议你关注一下wiki百科和SciPy，前 ...
refeihc 发表于 2010-5-13 02:06

只是供参考，原理都是非常简单的，具体实现就麻烦得多～
Intel的MKL就提供稀疏矩阵的BLAS库，但是还不完善。
wiki有时候会被撞墙的～
SciPy国内玩得人不多，目前为止Python科学计算相关的英文书籍不超过2位数，中文的话已经正式出版的还没有呢。与Matlab相比差很多。
又，Python只是胶水语言，纯python代码科学计算的效率惨不忍睹。。。把C/C++、Fortran代码打包封装，供Python调用才是上策。

pasuka

高效的求解器一直是有限元爱好者梦寐以求的东西，不过以本人的经验，
当前基于pc机的线性方程求解程序包还没有特别令人满意的。
本人遇到的问题和求解器现状：
1. 使用行压缩或列压缩存储刚度矩阵是必需的，对于超 ...
bbssbb 发表于 2010-5-13 06:56

使用PARDISO也觉得速度慢？
可以去水木社区的数值计算版看看，前几天刚讨论过这个问题。

s0101017

嗯 受益不少哇

refeihc

回bbssbb：

高效的求解器一直是有限元爱好者梦寐以求的东西，不过以本人的经验，
当前基于pc机的线性方程求解程序包还没有特别令人满意的。
......
我所期待的有限元pc机求解器：$ F! |; A! s1 I1 f6 K5 ~
1。具有通用性，可解对称，非对称，半正定，受条件数的影响较小。& Q) C8 r* ~/ I( A4 a
2。具有迭代法的效率。* K$ J, g% G% J( @  X
如果有这样的求解器存在，请各路数值线性代数的高手一定告知。
bbssbb 发表于 2010-5-13 06:56

我也关心高效的求解器，建议不时访问下列网站：
http://www.sccas.cn/gb/index.html
http://www.cs.sandia.gov/CRF/aztec1.html
http://www.cs.sandia.gov/
http://www.mcs.anl.gov/index.php
http://trilinos.sandia.gov/
http://www.mcs.anl.gov/petsc/petsc-as/index.html

这些机构就是吃这碗饭的。

refeihc

只是供参考，原理都是非常简单的，具体实现就麻烦得多～
Intel的MKL就提供稀疏矩阵的BLAS库，但是还不完善。
wiki有时候会被撞墙的～
SciPy国内玩得人不多，目前为止Python科学计算相关的英文书籍不超过2位数， ...
pasuka 发表于 2010-5-13 10:12

回pasuka：

同意你的观点，Python不适合大型计算，做script还是不错的。

refeihc

五、开源后的新较量：稀疏存贮格式与迭代算法

“三等车里拥挤得不堪了，头等车里只坐着三个人。中间不过隔了一扇玻璃门。愚蠢的人类呀，你们为什么不把这扇门打破了，大家坐得舒服些?”

郑振铎的这首“在电车上”形象生动地描述了生活中的一幅场景，你挤进了三等车箱，路上还有艰难行走的人群。你是该庆幸呢，还是该进一步向往舒服的头等车箱？结局取决于你能否打破一扇玻璃门。

我们不妨把开源软件比着三等车箱，首先想想，我们怎么能坐了进来。也就是问那么努力试图保持科技领先优势的美国人，怎么会允许我等轻松获得开源软件包呢？我觉得有几个原因：
(1) 一些开源软件在算法上已不具备先进性，别的国家开发出性能差不多的数学软件包难度不大，甚至已经开发出来了，所以闭源的意义不大，反而开源后有助于保住这块阵地。
(2) 世界上出现了一批反对软件垄断的潮流，让美国人觉得开源对自己也有一些好处，不仅可以获得好名声，还能扶持新生的力量。
(3) 有的开源软件所具有的传染性迫使其使用者不得不公开自己的源代码，这是一种隐形的制约，说明天底下没有免费的午餐。
(4) 希望我们安心地享用三等车箱——开源软件，他们呆在头等车箱可以舒服些。

这几方面我们很容易理解，但不会所有的人都乖乖地听话，于是更激烈的后续较量就展开了，那就说说这个话题吧。

涉及到矩阵求解器，首先要提到存贮格式和算法。矩阵的存贮格式有一般、等带宽、一维变带宽和稀疏格式等不同类型，每一类格式还可根据数据排列特点进一步细分，这里不详细讨论，不妨称前3种为非稀疏格式。矩阵算法分为2类，一类是直接算法，另一类是迭代算法。

存贮格式和算法可以有4种不同的搭配组合，尤其是“非稀疏存贮格式+直接算法”组合，目前被多数开源软件包采用。非稀疏存贮格式的主要不足之处是过于浪费存贮空间，以及太多不必要的计算，这限制了它被用于大型工程问题。

理论上，经典的“非稀疏存贮格式+直接算法”组合已趋于定型，这表明与此有关的开源软件包也不会有太多发展了。这并不是说就不能改进了，而是费力改进的油水不大，一方面在存贮空间上，如果使用等带宽存贮格式，不仅浪费空间，而且对前处理提出了一个至今无法彻底满足的要求——带宽优化；而如果使用一维变带宽存贮格式，虽然空间能省一些，但编程的麻烦程度一下增加了很多，而且带宽优化的要求仍然存在。何况不管采取哪种存贮方式，稀疏矩阵存贮方式所占用的存贮空间都要远远小于直接算法，请注意这里所说的“远远小于”是指：“随着问题规模越来越大，稀疏矩阵所占用的存贮空间数与直接算法所占的存贮空间数之比趋于零。”

在这样的背景下，“稀疏存贮格式+迭代算法”组合受到人们的广泛关注，目前这一方面的研究正在发展中。如果采用这种组合方式，计算中的许多困难都能得到克服，如
(1) 计算量大大减少，因为只操作非零元素。
(2) 后处理、杀死单元、变结构等方面的计算将变得方便，这是因为总刚度矩阵可以采用LIL格式存贮，实质上是没有把所有单刚矩阵元素加到总刚，而是分别保存。
(3) 约束条件的处理变得容易实现，对角元素改1法能轻松的去掉行和列，直接代入法也不再需要进行大数组搬家。
(4) 需要解代数方程时，再转换总刚成CSC、CSR或其它格式，计算量不大，存贮空间消耗量在数量级上也远低于矩阵直接解法的消耗量。
(5) 便于大型问题的并行计算。

不过稀疏矩阵存贮格式能否方便使用取决于能否实现高效率的迭代算法，目前世界上正对这一课题加紧进行研究，后续较量就在这一层面上展开，像我前面提到过的UMFPack就曾经被Matlab采用过，不过它协议最友好的时候也就是LGPL，后来又变成了GPL，这么变来变去，个中原因令人深思。

目前能够处理稀疏矩阵的开源软件有ATALAS、UMFPack、SuperLU、PETSc、MUMPS等，很难说谁就一定比其它的好。我当初选择UMFPack和SuperLU的主要原因是一种从众心理，相信别人的眼光，自己筛选起来挺劳神的。不选ATALAS是因为它是更早一代的代码。相信PETSc和MUMPS也一定很好，但在存贮格式和调用接口上它们似乎自成一派，看不清未来的走向，因此没有选它们。

关于大型稀疏矩阵的迭代计算方法，近期已有许多成果，最新的如Arnoldi、PCG、Krylov、GMRES方法等，我没有钻研过。以前学习过的方法有幂、反幂、松驰、最速下降等，知道迭代法如果能快速收敛，效率是非常惊人的。而如果收敛慢，甚至于是不收敛，那就只能望着干着急。

常言道“三十年河东，三十年河西”，激烈竞争的软件行业可能风云变化更快，没准河东变到河西三年都不用。开发软件的人总是希望用最好的开源软件包，但又说不准今天的最好明天会落伍。所以就希望能有一个一劳永逸的办法，万一有朝一日有后来居上者，我只须把接口一换，其它代码不用重写，岂不是好。可惜的是上述开源软件开发者都自行其是，各干各的，你如果用了其中一种，想换另一种，那对不起，有关矩阵生成、程序调用的代码请你重写。

于是，就提出了另一个要求，选择矩阵求解器时，还要考虑其存贮格式和调用方式的通用性，这是求解器开发要关心的又一个关键问题。

caoer

GMRES 挺好，preconditioner功能非常强劲

refeihc

楼上为什么感兴趣PETSc呢？

其实，我也挺感兴趣数值算法，但是时间、精力、能力都有限，对软件包的一点点尝试，已经感觉自己渺小了

caoer

个人的偏好把，嗬嗬，主要是他的并行做的很好，最近推广的也很好

hillyuan

29# refeihc

: 目前的开源软件包多数采用直接算法，少部分基于稀疏矩阵存贮方式和并行迭代计算的开源软件包还在发展中
=> Are you sure?
   Maybe we cannot wish a great development of direct solver, at least in near furture, because direct solver can solve EVERYTHING theoretically. However, iterative may fail to converge or too slow to converge. A typical example in FEM analysis is when shell element exists, we generally don't use iterative solver. Researchers need much to do to improve the effiiency of iterative solver, generally by adopt a good precondition calculation.
  Some new researches concerns parrellel calculation, include GPU parrellel calcuation of iterative solvers.

  You need clearify the "sparse matrix" and "iterative solver" (the so-called band matrix in nonsense in application now, both for direction and iterative solvers) used in your post. For examples, UMPack is a direction, not iterative, solver using sparse matrix. All the solver list in your post adopt sparse matrix. Just forget band matrix from now.

refeihc

33# hillyuan

...
涉及到矩阵求解器，首先要提到存贮格式和算法。矩阵的存贮格式有一般、等带宽、一维变带宽和稀疏格式等不同类型，每一类格式还可根据数据排列特点进一步细分，这里不详细讨论，不妨称前3种为非稀疏格式。矩阵算法分为2类，一类是直接算法，另一类是迭代算法。

存贮格式和算法可以有4种不同的搭配组合，尤其是“非稀疏存贮格式+直接算法”组合，目前被多数开源软件包采用。非稀疏存贮格式的主要不足之处是过于浪费存贮空间，以及太多不必要的计算，这限制了它被用于大型工程问题。

理论上，经典的“非稀疏存贮格式+直接算法”组合已趋于定型，
...
refeihc 发表于 2010-5-16 20:16

楼上说得对，应该有4种组合情况，以前表述没有加以区分，是个失误，现在修改了；UMFPack就不补充介绍了。请再看看有什么不妥。

谢谢你指出问题！

pasuka

五、开源后的新较量：稀疏存贮格式与迭代算法

“三等车里拥挤得不堪了，头等车里只坐着三个人。中间不过隔了一扇玻璃门。愚蠢的人类呀，你们为什么不把这扇门打破了，大家坐得舒服些?”

郑振铎的这首“在电车上 ...
refeihc 发表于 2010-5-16 20:16

“稀疏存储格式+迭代算法”组合受到关注？
首先，这个表述有点拗口。直接说：大规模稀疏矩阵的迭代求解器
关注肯定是有不少，主流商业CAE软件目前还是直接求解器和迭代求解器并存，且直接求解器更成熟实用些。
其次，2维等带宽和1维变带宽实际在大规模科学计算中已经予以扬弃，虽然国内很多有限元教科书还是乐此不疲的在介绍。。。
再者，矩阵重排漏讲了，这已经不是简单的半带宽最小那么简单，AMD、RCM等等算法（具体这些算法的简单介绍，参考Matlab的symamd、symrcm函数的帮助文档，主要作用，用一句话概括：矩阵分解过程中，零元位置不出现非零元），任何直接求解器必然会有矩阵重排优化这一步骤，通常会调用例如METis这样的开源软件包。
btw，更关注稀疏矩阵广义特征值的求解，lz提到的幂法、反幂法都是教科书上介绍特征值计算的方法吧？

refeihc

呵呵，咱要不抛块土砖，哪能引来恁多宝玉。

我写这个帖子感觉有些吃力了，辛苦倒在其次，主要是知识储备明显不够。但是各位高手加入进来讨论，带给我意外的收获。所谓：朝闻道，夕死可矣！所以后面还得咬紧牙关坚持。

楼上，谢谢补充！
遗漏是无法避免的，我们现在能看到的有限元、数值算法等教科书既有国内的，也有国外的，我印象中，国内的这些书基本不会详细介绍大型稀疏矩阵迭代算法，反而介绍幂法等基础算法较多。

如果哪位高手了解到目前国外在这方面出了比较牛的教材，涉及大型稀疏矩阵方程迭代求解、迭代求广义特征值、求广义复特征值等内容，希望能在这里推荐。

refeihc

六、ublas——统一数据格式？
易经语：易穷则变，变则通，通则久。我们的老祖宗非常重视变通之理，认为只有变通才能长久。当别的开源软件开发者都各行其路时，唯有变通才能立于不败之地。

前面说过选择矩阵求解器时，不仅要注意它的当前性能，还要考虑其存贮格式和调用方式的通用性，这其实是一种变通的思想。对于数学软件包的开发者来说，重视这一问题有助于所开发软件的长期健康发展，而对于使用者来说重视这一问题可以少走弯路，获得源源不断的技术支撑。实际上，这已经不是简单地只关心软件包的现状了，还要重视其未来的发展趋势，目前网上的讨论较少涉及到这一方面。

一个软件包的未来发展如何，是兴盛还是衰落，大多数人是不会关心的。但对于那些希望站在一个较高平台上进行软件开发的人来说，则不能不重视。

如果所选用的软件包一直傲视群伦，技术上保持领先优势，那应该庆幸当初的眼力不错。这使得我们可以只关心自己的软件功能，别的什么都不用操心，开源软件包就这么一直用下去，可谓省心省力。

但软件行业风云变幻，谁也料不定将来会怎么样，万一哪天江湖座次重新排过——这很可能就是一年半载的事，原来所用的软件包落伍了，更好的开源软件包出现，此时当然想更换一下软件包。这时得有思想准备，必须重新编写调用新软件包的有关代码，通常这样的修改还能接受，毕竟工作量不算大。但某些情况下，可能连数据结构都要重新设计，这可能导致程序许多模块都要重新编写，这样的情况是开发者最不愿意看到的。

怎么办？我们可能预见不到将来，但是我们能够尽量规避未来的不利情况。这里介绍一下我了解到的又一个开源软件包ublas，它对于上面的问题比较重视，目前正在发展中。

熟悉C++的人都比较了解boost，这里不打算专门介绍它。虽然boost中还有许多的其它数学软件包（如math/special_functions、math/statistical distributions等），但与矩阵计算有关的数学软件包只有ublas。国内可能讨论其它软件包(如lapack等)的人更多，但ublas还是有些特殊的地方值得关注。

1 首先，boost如同它的名字一样带给人强劲有力的印象。人们普遍认为它是继STL之后C++委员会陆续完成的一系列大动作。矩阵类以及相关的数学计算一直受到业界的广泛关注，人们确信Fortran更擅长大量的科学计算工作，而面向对象的C++语言则在体系架构方面具有显著优势，因此做为boost家庭的一员，ublas产生的目的就是要兼顾这两种语言的优点。

2 其次，ublas为使用者提供了最基本的矩阵数值计算（类似于blas的三个水平），这可满足大多数情况下的要求。

3 对于更高的计算要求，ublas可通过boost.numerical.bindings来拓展功能，从而建立起一座通向多个不同数学软件包（目前有amos、atlas、blas、lapack、mumps和umfpack）的桥梁。

4 最后，它似乎目标很大，别的软件包只想树起一面旗帜，而它可能更想一统江湖。我这样说是因为它的矩阵类实现了基于算子重载的数学自然语言语法，而不同属性矩阵的交叉计算则是采用了基于矩阵表达式模板(expression template)的调用方式，这样一来其矩阵的数据存贮方式达到了一定程度的通用性。

先看一个例子。

例1  MTL要通过下面的语句
  template <class T, class Shape=rectangle<>, class Storage=dense<>, class Orientation=row_major>
  struct matrix {
    enum { Shape_id = Shape::id, Shape_M = Shape::M, Shape_N = Shape::N };
    //: The generated type
    typedef typename IF< EQUAL< Shape_id, RECT>::RET,
                        typename generate_rect<T,Orientation,Storage,
                                               Shape_M, Shape_N>::RET,
                    typename IF< EQUAL< Shape_id, DIAG>::RET,
                        typename generate_diagonal<T,Storage,
                                               Shape_M, Shape_N>::RET,
                    typename IF< EQUAL< Shape_id, BAND>::RET,
                        typename generate_banded<T,Orientation,Storage,
                                               Shape_M, Shape_N>::RET,
                    typename IF< EQUAL< Shape_id, TRI>::RET,
                        typename generate_triangle<T,Shape,
                                               Orientation,Storage,
                                               Shape_M, Shape_N>::RET,
                    typename IF< EQUAL< Shape_id, SYMM>::RET,
                        typename generate_symmetric<T,Shape,
                                              Orientation,Storage,
                                               Shape_M, Shape_N>::RET,
                        generators_error
                        >::RET
                        >::RET
                        >::RET
                        >::RET
                        >::RET type;
  };
来声明一个模板化的Matrix结构，我相信任何一个人当看到这样麻烦的矩阵类型时，虽然也会佩服程序设计者将typedef和宏语句用得出神入化，但会感觉到这种数据结构只能用于MTL自己，如果想换其它的软件包，势必有大量的代码需要重写，因此可以预料MTL后面的路不会太好走，除非它的开发队伍像苹果公司那样，能我行我素地闯出一条独具特色的路。

而ublas的作法与MTL完全不同，从下面的2个例子可以看出它打算基于统一的数据格式来实现对其它数学软件包的调用。

例2  调用lapack的方式
    ublas::matrix<double, ublas::column_major> a(5,5),b(5,2);
    ......  // 生成矩阵数据
    lapack::gesv(a,ipiv,b);  // 调用lapack求解方程

例3 调用UMFPack
    ublas::compressed_matrix<double,ublas::column_major,0,ublas::unbounded_array<int>,ublas::unbounded_array<double>> A(5,5,12);
    ublas::vector<double> B(5),X(5);
    ......  // 生成矩阵数据
    umf::symbolic_type<double> Symbolic;
    umf::numeric_type<double> Numeric;
    umf::symbolic (A,Symbolic);
    umf::numeric (A,Symbolic,Numeric);
    umf::solve (A,X,B,Numeric);

可以看出，不管调用什么软件包，数据格式都是ublas的，在这方面它应该做了大量的研究。ublas今后能走多远，我们拭目以待。

Mesh_Li

写的真好 学习了 谨受教！

pqjlmy

今天一口气读完了这楼的所有帖子，唯一感受就是自己的无知，自己知识的可怜，希望楼主多发些这里帖子！同时顺便问一句：《矩阵计算/中国科学院研究生教学丛书》作者: （美）G.H.戈卢布（Golub，G.H.），（美）范洛恩（Van Loan，C.F.） 著，袁亚湘 等译
这本书是不是对学习lapack有很大帮助？？？

refeihc

楼上的兄弟，袁亚湘的这本书我只看过目录，我当初学习Lapack是直接看的帮助文件，所以你的问题我回答不了，不过我相信会有高手回答的。