GEO_WU

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日星期一

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

2013年4月18日星期四

latex宏包安装

在ubuntu下安装latex宏包过程如下：

下载宏包.sty文件或者是自己编写的.sty文件；
找到系统默认宏包的位置，通常是在/usr/share/texmf/tex/latex目录下；
将下载的宏包或自己编写的宏包移到系统对应的宏包位置；
更新一下： sudo texhash

举例：

自己编写的beamercolorthemeLUH.sty文件；
系统对应的beamercolortheme宏包位置为：/usr/share/texmf/tex/latex/beamer/base/themes/color
sudo mv beamercolorthemeLUH.sty /usr/share/texmf/tex/latex/beamer/base/themes/color
sudo texhash

2013年4月17日星期三

关于字符编码的一点知识

分享两篇介绍字符集和字符编码的博文：

此外，常用的latin字符集介绍见另一个链接：

https://zh.wikipedia.org/wiki/ISO/IEC_8859

2013年4月4日星期四

英文科技文献写作中的一些细节

在英文科技文献写作过程中，经常在一些细节上模棱两可，这里将这些细节做一汇总。

1. 参考文献和引用的撰写（基于journal paper）
不同的期刊有不同的参考文献格式，并没有一个统一的参考文献格式，这里只给出一种常用的格式。参考文献包含如下信息：作者+年份+标题+期刊名+卷号+页码，这里重点讨论作者的写法。首先列出常用的几条准则：

所有参考文献是按第一作者的last name (family name, or surname)的首字母顺序排列；
姓名缩写的规则是last name （用于正式场合）不缩写，given name （first name and middle name）可以缩写；
参考文献中必须包含所有在文章中出现的引用，and no more!
同一作者在一年里有多篇文章，需要在年份后加小写字母来区分这些文章；
直接引用别人的文章，必须加引用（但不推荐这种引用方法）；转述别人文章中的idea或者facts，也必须加引用（大多数为这种情况）。
不仅仅book和paper, internet source, computer software, e-mail and verbal communication都需要加以引用说明;
用了别人文章中的图，也必须加以引用；

接下来给出常用的参考文献和其对应的引用：

Single Author Bugjuice, B. 1970. Physiological ... (Bugjuice, 1970)

Two Authors Timm, T., Bugjuice, B. 1990. The role ... (Timm and Bugjuice, 1990)

Muitiple Authors Bugjuice, B., Timm, T., Cratchet, R. 1990. The ... (Bugjuice et at., 1990)

另外两种常见情况：

同时引用多个文章，如（Gumward, 1952; Bugjuice, 1970）。一个作者一年有多篇文章(Bugjuice, 1970a)

注：在写引用时，只给出last name就可以了。

参考资源：

http://abacus.bates.edu/~ganderso/biology/resources/writing/HTWcitations.html

http://tim.thorpeallen.net/Courses/Reference/Citations.html

2. 数字和单位之间是否要加空格？ Yes.

3. 单位符号是用斜体还是正体表示？正体。

4. 括号前后是否需要加空格？ Yes. 如： We spent NT$3,000 (about 100 U.S. dollars) on dinner. 但和括号里的第一个和最后一个字之间不加空格。

5. i.e., e.g. 和etc.的用法
i.e. 是拉丁文 id est 的缩写，它的意思就是“那就是说，换句话说”，等同于“that is，in other words” ，目的是用来进一步解释前面所说的观点。
e.g. 是拉丁文 exempli gratia 的缩写，它的意思是“举个例子，比如”，等同与“for example”,目的就是用几个例子来说明前面的观点。
etc.就比较好理解了，它是 etcetera 的缩写，意思是“等等”，相当于“and so on”
e.g. 和 etc. 不能出现在同一句话中，因为 e.g. 是表示泛泛的举几个例子，并没有囊括所有的实例，其中就已经包含“等等”，如果再加一个 etc. 就画蛇添足了，例如下面这句话就是错的：
Writing instructors focus on a number of complex skills that require extensive practice (e.g., organization, clear expression, logical thinking, etc.)
注意，用它们的时候，不要把“小点”给落了~

参考资源：http://bbs.ednchina.com/BLOG_ARTICLE_43177.HTM

2013年3月25日星期一

ubuntu下安装Lapack

1. 关于Lapack和Blas
Lapack (Linear Algebra Package)是一个以Fortran90编写，用于数值计算的函数集。例如，用于解线性方程组、计算特征向量、计算矩阵QR分解、奇异值分解等问题。
Blas(Basic Linear Algebra Subprograms)，基础线性代数程序集，用于基本的矢量和矩阵运算。最初发布于1979年，并用于建立更大的数值程序包（Lapack）。一级Blas负责矢量-矢量运算，二级Blas负责矢量-矩阵运算，而三级Blas负责矩阵-矩阵运算。由于Blas的高效，在高质量线性代数软件中常常使用它。

2. Lapack安装
1> 下载Lapack
http://www.netlib.org/lapack/ 最新版本3.4.2
2> 解压
tar -zxf lapack-3.4.2.tgz
3> 编辑make.inc文件
进入解压后的文件夹，会有一个make.inc.example文件，将该文件重命名为make.inc文件；并查看该文件，默认编译器为gfortran，根据需要修改对应的参数。在ubuntu下，采用gfortran作为编译器，一般不许要修改该文件。
4> 编辑Makefile文件
将下列语句（设定生成的库，默认不包含blaslib）
lib: lapacklib tmglib
# lib: blaslib variants lapacklib tmglib
改成：
# lib: lapacklib tmglib
lib: blaslib variants lapacklib tmglib
生成的库中包含了blaslib。
5> make
安装并测试Lapack。
6> 将生成的静态库(.a)移到/usr/lib文件夹下
生成的静态库包含libblas.a liblapack.a libtmglib等三个库。
sudo mv *.a /usr/lib

3. Lapack使用
在编译程序时，语句如下（其中test.f90文件中包含了lapack里的routines）：
gfortran -o test test.f90 -llapack -lblas

4. Lapack概述
1> Levels of Routines
Lapack的程序分三个层次，分别是driver， computational和auxiliary。driver routines解决一个具体的问题，如求解线性方程组或计算一个实对称阵的特征值。如果用户的问题刚好可以用该层次的程序解决，推荐用该层次的程序。这类程序通常调用几个computational routines或者一些其他的子程序。computational routines解决单独的一个计算任务，如LU分解。用户可以根据自身的需求调用该类程序。auxiliary routines这里不做阐述。一般用户一般调用driver和computational routines较多。
2> 数据类型和函数命名规则
Lapack函数命名格式为XYYZZZ，其中：
X表示数据类型：

S 表示单精度；
D 表示双精度；
C 表示复数；
Z 表示双精度复数；

YY表示矩阵类型：

DI 表示对角阵；
GE 表示一般矩阵；
GT 表示一般三角阵；
OR 表示正交阵；
PO 表示正定阵；
TR表示三角阵；

ZZZ表示具体的计算问题：

SV或SVX（求解线性方程组），前者是simple版，后者是expert版（输出一些中间结果，耗费内存更多）
QR（QR分解）
LU （LU分解）
...

3> 几个常用的函数名

dgesv 求解线性方程组AX=B
dgels 线性方程组AX=B的LS解
xyyTRF 将矩阵yy进行三角分解（根据矩阵yy的类型，选择不同的分解方法，如GE矩阵用LU分解，对称正定矩阵采用Cholesky分解）
xyyTRS 基于xyyTRF分解的结果求解AX=B
xyyTRI 基于xyyTRF分解的结果求解$A^{-1}$

因此，dgesv等价与dgetrf和dgetrs的组合，dposv等价与dpotrf和dpotrs的组合，dgels等价与dgeqrf和dtrtrs的组合。

更多内容参见Lapack User's Guide, http://www.netlib.org/lapack/lug/

5. 不错的资源

官方手册 http://www.netlib.org/lapack/lug/
用法举例 http://www.nag.co.uk/lapack-ex/lapack-ex.html#begin
中文帮助手册 http://wenku.baidu.com/view/d6c95d333968011ca30091fc.html
程序查询 http://www.netlib.org/lapack/explore-html/index.html

补充：
Lapack是用fortran写的，Lapacke是它的C语言接口，先安装Lapack,然后安装Lapacke,那么在程序中可以直接调用C函数来实现需要的功能。

订阅：博文 (Atom)