GEO_WU: 五月 2013

2013年5月30日星期四

自定义colormap meep (zt)

matlab 自定义colormap meep (from red through white to blue)

I have been thinking of using the colormap used in the MEEP software (an open source FDTD code) as a nice alternative to the classical Matlab colormap. Here is the short Matlab code to generate it and the resulting images:


mincolor    = [1 0 0]; % red
mediancolor = [1 1 1]; % white   
maxcolor    = [0 0 1]; % blue      

ColorMapSize = 4;
int1 = zeros(ColorMapSize,3); 
int2 = zeros(ColorMapSize,3);
for k=1:3
    int1(:,k) = linspace(mincolor(k), mediancolor(k), ColorMapSize);
    int2(:,k) = linspace(mediancolor(k), maxcolor(k), ColorMapSize);
end
meep = [int1(1:end-1,:); int2];

colormap(meep); 
colorbar

参考资料:

http://meyavuz.blogspot.de/2011/02/meep-colormap-for-matlab.html

2013年5月9日星期四

线性方程组求解与矩阵求逆

1. 线性方程组
在大地测量数据处理方面，线性方程组的求解无疑是最重要的问题之一。线性方程组的通用形式可写为

$AX=B$ ，系数矩阵A的维数为

$(M，N)$ ，

$M$ 为观测方程的个数（等于观测值的个数，1个观测值一个方程），

$N$ 为需要求解的系数个数，因此向量

$X$ 的长度为

$N$ ，观测值向量

$B$ 的长度为

$M$ 。根据

$M$ 和

$N$ 的大小观测，可以将线性方程组分为三类：

$M > N$ , 超定方程组（也称矛盾方程组），此时通常求方程的LS解；
$M = N$ ，适定方程组，有唯一解；
$M < N$ ，欠定方程组，有无数组解，但通常求其最小范数解。

接下来主要讨论超定方程组和适定方程组的情况。

对于超定方程组，根据LS准则，可得到法方程：

$A^TAX=A^TB$ ，令

$C=A^TA，D=A^TB$ ，则

$C$ 的维数为

$(N，N)$ ，D的维数为

$(N，1)$ ，此时法方程的求解与适定方程组一致。下面的讨论以适定方程组为例。

2. 线性方程组的求解与求逆的关系

线性方程组

$AX=B$ 的求解在数学上表示为

$X=A^{-1} \times B$ ，但是在实际数值计算过程中，并不是一定要去计算

$A^{-1}$ ，如常用的消去法，三角分解前推回代法等。

矩阵求逆等价于求解方程组

$AX=I$ ，理论上与求解线性方程组一致。

在数值计算过程中，我们必须考虑计算的速度以及对内存的耗用情况，此时求逆和线性方程组的求解是两个不同的概念。

在求解线性方程组之前，必须明确你所需要的结果。如果只要最后结果

$X$ ，则尽量不要去求逆；

如果系数阵

$A$ 的求逆无法避免，尽量选择快速的求逆算法。当然，对于小矩阵，就不用考虑两者的差别了。

用matlab函数表示以上情况：

$X=A\B$ (基于LU分解，前推回代直接求解参数向量

$X$ )；

$X=A^{-1} \cdot B$ （基于求逆求解参数向量

$X$ ）；

通常第一种方法速度快，占用内存少，但无法给出

$A^{-1}$ 的信息。

3. 线性方程组的常用数值解法

通常把线性方程组的数值解法分为直接法和间接法。所谓直接法，就是经有限次运算即可求得方程组的精确解的方法；迭代法将求解的方程组转化为一个无限序列，其极限就是方程组的解。直接法一般来说用于求解系数矩阵阶数不是太高的问题，工作量较小，精度较高，但程序较复杂；迭代法主要用于一些高阶问题，一般程序较为简单，但工作量有时较大，耗用内存较多。

1> 直接法

a. 克莱姆法则

$x_i=\frac{det(A_i)}{detA}$

算法复杂度：

$(n+1) \times [(n-1) \times n! + n!]$

运算量非常大，理论上可行，实际计算中很少用；

三角形方程组的求解简单快速，所以通常将系数阵A转化为三角阵再进行运算。不同的分解形式对应不同的分解方法。求解下三角阵L的过程，称为前推

$(X_1-X_2---X_n)$ ;求解上三角阵U的过程，称为回代

$(X_n---X_2-X_1)$ 。前推和回代的计算量一样，为

$\frac{n^2}{2}$ 次乘法和加减法以及n次除法。

b. 高斯消去法

利用消去法，将方程组变换为上三角形方程组；整个消去过程需完成

$O(\frac{1}{3}n^3)$ 次乘法和加法运算；回代过程需

$\frac{n^2}{2}$ 次乘法和加减法以及n次除法，故用消去法求解n阶方程组的计算量约为

$O(\frac{2}{3}n^3+2n^2)$ ；

从高斯消去法进一步完善，产生了列主元和全主元消去法，这两种方法在消去过程中的计算量进一步增加，但方法上更严谨。

每一步消去过程对应系数阵A左乘了一个初等变换矩阵，这一系列初等变换矩阵的乘积对应一个下三角矩阵L，故高斯消去法理论上等价与LU分解法。

c. LU分解法

将系数阵A直接分解成两个三角阵的乘积，根据矩阵运算的规则计算每个三角阵的元素。分解后

$AX=B$ 的计算转化为两个三角形方程组的求解：

$LY=B$ （前推，计算出Y）

$UX=Y$ （回代，计算出X）

算法复杂度为

$O(\frac{2}{3}n^3+2n^2)$ ，所以直接分解法与消去法的计算工作量基本上是相同的。

求解逆阵的算法复杂度为

$O(2n^3)$ ，是求解线性方程组的3倍。

d. 正交分解法

将矩阵A分解为一个正交阵和一个上三角阵的乘积，

$A=QR$ 。根据分解方法，通常有Householder方法（镜像映射法）和Schmidt正交化方法。

运算复杂度为：

$O(\frac{4}{3}n^3+2n^2)$ ，与上面两种方法相比，运算量多一倍。

该方法在求解超定方程组时，较为优越，原因如下：

$A^TAX=A^TB$ with

$A=QR$ , then

$R^TQ^TQRX=R^TQ^TB$ and finally,

$RX=Q^TB$ ，因此，QR分解大大减少了计算量。

e. 对称正定矩阵的平方根法和 $LDL^T$ 分解法
对称正定矩阵是特殊矩阵，可以用平方根法和

$LDL^T$ 分解法（通常也叫做Cholesky decomposition），其运算量和存储量较普通消去法和LU分解法均节省一半，为

$O(\frac{1}{3}n^3+2n^2)$ 。
平方根分解：

$A=LL^T$

$LDL^T$ 分解：

$L$ 为单位下三角阵

2> 迭代法
对于阶数不高的线性方程组，直接法有效，但对于阶数很高的矩阵，有时需要迭代法来解决。迭代的规则不同，迭代法也就不同。通常分为：线性迭代和非线性迭代；一阶迭代和P阶迭代；定常迭代和非定常迭代。具体，可参考教科书，如“数值计算方法（冯康）”。
共轭梯度法：非线性迭代法，适用于对称正定阵，不需选取任何迭代参数，使用方便。
此外，对于高阶方程组的处理通常还采用分块矩阵处理方法。节省内存需要量，但增加了计算量，需要更多的计算时间（由于内外存交换数据需要时间）。

4. 矩阵求逆
除了必须要求矩阵的逆，否则尽量不采用求逆阵的方法来求解线性方程组。
求逆阵的算法复杂度为

$O(2n^3)$ ，大约是同等方法求线性方程组

$O(\frac{2}{3}n^3)$ 的3倍。
矩阵乘积的算法复杂度和矩阵求逆是等价的。
矩阵求逆当前最快方法的复杂度为

$O(n^{2.37})$
高斯消去法，LU分解法，SVD等，复杂度都为

$O(n^3)$ ，差别在于前面的系数。

特殊矩阵的处理通常会有特殊的方法，能大大减少运算量和存储量。

2013年5月7日星期二

lapack 线性方程组最小二乘求解函数dgels

1. DGELS函数简介

DGELS函数是利用最小二乘求解线性方程组（超定或欠定）AX=B或A'X=B的函数，处于driver level。利用LU分解或者QR分解进行求解。

2. DGELS函数调用格式

DGELS(TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK, INFO)

输入参数：

TRANS（CHAR）：值为'N'，表示No transpose, 对应AX=B；值为'T',表示Transpose，对应A'X=B;

M,N（INT）: 系数矩阵A的维数，M>N,超定（下面只以超定方程组为例）,M<N，欠定；

NRHS（INT）：观测值向量B的列数（通常为1）；

A（Double）：系数阵；

LDA（INT）：同DGESV函数中的LDA，LDA>=max(1,M)，通常等于M；

B（Double）：观测值向量，当TRANS=‘N’时，B（M，NRHS），当TRANS=‘T’时，B（N，NRHS）；

LDB（INT）：同DGESV函数中的LDB，LDB>=max(1,M)，通常等于M；

LWORK（INT）：Length of Work, 数组WORK的长度，见http://www.netlib.org/lapack/lug/node118.html；

输出参数：

B（Double）：B(1:N)为计算的结果；B（N+1：M）误差信息，具体见函数说明；

A（Double）：M>N, QR分解矩阵信息；M<N, LU分解矩阵信息；

WORK（Double）：维数为LWORK，存储空间(WORKspace)信息，具体见说明http://www.netlib.org/lapack/lug/node120.html

INFO：=0，成功执行；=-i，第i个参数illegal；=i，A的第i个对角元素为0，A非满秩阵，无LS解。

3. DGELS函数举例：

http://www.nag.co.uk/lapack-ex/node45.html#DGELS

4. 同类型函数

DGELSY：利用完全正交分解计算AX=B的LS解
DGELSS：利用奇异值分解计算AX=B的LS解
DGELSD：利用奇异值分解（divide and conquer）计算AX=B的LS解

5. 参考资料

2013年5月6日星期一

lapack 线性方程组求解函数dgesv介绍

1. DGESV函数简介
DGESV函数是lapack函数库求解线性方程组AX=B的函数，处于drive level。利用LU分解进行求解，A=PLU，P为置换矩阵，L为单位下三角矩阵，U为上三角矩阵。

2. dgesv函数参数介绍
a) 形式
DGESV(N, NRHS, A, LDA, IPIV, B, LDB, INFO)
输入参数：
N（INT）：线性方程组的数目，矩阵A的维数；
NRHS（INT）：矩阵B的列数（一般为1，单列）；
A（Double）：系数矩阵A，维数（LDA，N），通常LDA=N；
LDA（INT）：Leading Dimension of A, LDA>=max(1,N);
B（Double）：矩阵B，维数（LDB，NRHS），通常LDB=N，NRHS=1；
LDB（INT）： Leading Dimension of B，LDB>=max(1,N);
输出结果：
A（Double）：存储LU分解的L和U矩阵，L为单位下三角阵且对角线的单位值1没有存储；
IPIV（INT）：维数为N，利用该矩阵的值可以推求置换矩阵P，starting from a unit matrix E, E矩阵的第i行与E矩阵的第IPIV（i）行进行交换。
如IPIV=[2 2 3 4]，则P矩阵的结果为：

B（Double）：存储最终的求解结果X值；

INFO（INT）：求解状态参数。INFO=0，成功执行；INFO<0，则INFO=-i，第i个参数是illegal value；INFO>0, U(i,i)=0, LU分解成功执行，但是U是奇异的，最终结果无法计算；

DGESV函数调用了Lapack Computation level的两个函数DGETRF和DGETRS，前者是进行LU分解的函数，后者是利用分解的结果进行前向和后向推求最终结果X。简单介绍下这两个函数的调用格式：

DGETRF（M，N，A，LDA，IPIV，INFO），矩阵A的维度为（M，N），其它参数同上，返回结果为A（PLU），IPIV以及INFO。

DGETRS（TRANS，N，NRHS，A，LDA，IPIV，B，LDB，INFO），求解AX=B或A'X=B，该函数只有在INFO=0时才可以进行。输入变量'TRANS'为字符型，其意义如下：'N'对应AX=B； ‘T’对应A'X=B; 'C' 对应A'X=B (conjugate transpose)，其它变量的意义同DGESV函数。

3. dgesv函数举例
http://www.nag.co.uk/lapack-ex/node5.html
最终返回结果：
B=X= [ 1.0000 -1.0000 3.0000 -5.0000];
A=PLU=[
5.2500 -2.9500 -0.9500 -3.8000
0.3429 3.8914 2.3757 0.4129
0.3010 -0.4631 -1.5139 0.2948
-0.2114 -0.3299 0.0047 0.1314]

IPIV=[2 2 3 4]; 从一个单位阵开始，单位阵的第一行和第2行交换，第二行和第2行交换，第三行和第3行交换，第四行和第4行交换，

所以：

P=[0 1 0 0; 1 0 0 0; 0 0 1 0; 0 0 0 1];

L=[1 0 0 0; 0.3429 1 0 0; 0.3010 -0.4631 1 0; -0.2114 -.3299 0.0047 1];

U=[5.2500 -2.9500 -0.9500 -3.8000; 0 3.8914 2.3757 0.4129; 0 0 -1.5139 0.2948; 0 0 0 0.1314]

经过验证PLU=A；

4. dgesvx函数
dgesvx是expert dgesv函数，该函数可以返回更多的中间结果，如计算矩阵A的条件数，估计解的误差范围等等，相应的，函数也增加了很多参数，所以占用的内存也就更多。具体使用可以参见user book http://www.netlib.org/lapack/lug/node26.html和应用举例 http://www.nag.co.uk/lapack-ex/node6.html#DGESVX

2013年5月4日星期六

fortran编译器

Windows和Linux下Fortran编译器众多，这里简要介绍一些常用的编译器。

一. Windows平台下
在Win平台下，用的较多的是基于微软Visual Studio而集成的一整套IDE、商业编辑器：CVF和IVF。
1. CVF（Compaq Visual Fortran）
最常见的是CVF6.5，基于VS6，是win98时代的东西，更新到V6.6c被Intel收购了，即后来的IVF。CVF只支持到Fortran 95标准。
2. IVF（Intel Visual Fortran）
CVF的后继者，支持扩展的2003标准。IVF对Fortran77，90，95的兼容性不错。

二. Linux平台下
Linux下的fortran编译器很丰富，而且几乎都是免费的，如g77, gfortran, ifort 以及g95。
g77以前是GNU里的成员，后来被gfortran替代了。g95是在gfortan的发展过程中独立出去的，现在就只有作者一个人在维护。ifort是Intel的产品，对于单独的用户是免费的。这四个编译器用的较多的是gfortran和ifort，关于这两个编译器的讨论也集中在对老版fortran代码的兼容性（容错性）以及编译后的程序运行效率问题。总体来说，ifort的兼容性更好些，效率问题不太清晰。以下下面这两篇帖子对这个问题做过一些讨论，希望读完后自己能有一个总体的评判。

注：关于ifort在ubuntu下的安装可以参加另一篇博文：

http://geowu.blogspot.de/2012/07/ubuntuifort-xe.html

订阅：博文 (Atom)

Who am I